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Abstract 


Subsystem Density-Functional Theory (DFT) is an emerging tech¬ 
nique for calculating the electronic structure of complex molecular and 
condensed phase systems. In this topical review, we focus on some re¬ 
cent advances in this field related to the computation of condensed 
phase systems, their excited states, and the evaluation of many-body 
interactions between the subsystems. As subsystem DFT is in princi¬ 
ple an exact theory, any advance in this field can have a dual role. One 
is the possible applicability of a resulting method in practical calcula¬ 
tions. The other is the possibility of shedding light on some quantum- 
mechanical phenomenon which is more easily treated by subdividing a 
supersystem into subsystems. An example of the latter is many-body 
interactions. In the discussion, we present some recent work from our 
research group as well as some new results, casting them in the current 
state-of-the-art in this review as comprehensively as possible. 
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1 Introduction 


1.1 The idea behind subsystem DFT 

The founding idea of subsystem Density-Functional Theory (DFT) is that 
the electron density of a system can be represented by the sum of the electron 
densities of a collection of subsystems. Namely, 

Ns 

P(r) = 5^P/(r), (1) 

I 

where Ns is the number of subsystems chosen. To ensure iV-representability, 
each subsystem electron density is constructed to integrate to a predeter¬ 
mined and not necessarily integer number of electrons, Nj, such that their 
sum recovers the total number of electrons, e.g. '^jNj = Ne- 

Our research group, as well as others, often mention the above equation 
as a means to simplify the electronic problem by invoking the principle of 
divide-and-conquer. In fact, in the absence of underlying approximations, the 
above equation is hardly useful. One would be inclined to think that Eq. (1) 
conveys the idea that a system (supersystem, hereafter) can be calculated 
more efficiently when the problem is fragmented into smaller sub-problems. 
However, the subsystem electron densities on the r.h.s. of Eq. (1) are generally 
unknown, and so is the supersystem density, p. Thus, from the point of view 
of quantum mechanics, a density partitioning should not help. 

An important realization is that the density partitioning in Eq. (1) is 
not unique in the absence of additional requirements (constraints). It is 
the nature of those constraints that will make it possible to apply sensible 
approximations which can result in an efficient computation. For example, if 
the subsystem electron densities can be constructed such that their overlap 
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is minimized, algorithms able to speed up the computation by exploiting 
locality of the electronic structure [1] can be applied. 

There is another approach which is philosophically orthogonal to what 
is mentioned above - it prescribes foregoing the quest to computational ef¬ 
ficiency, and pursuing another goal. Namely, hnding the most chemically 
meaningful partition of the supersystem into subsystems. Since the pioneer¬ 
ing works of Bader [2,3], Hirshfeld [4], and Parr et al. [5,6] this quest has 
drawn the attention of many research groups [7-10]. Such interest stems 
from the historical fact that chemists gain chemical insight from properties 
of atoms (or fragments). Thus, a central problem of chemistry should be hnd¬ 
ing an appropriate way of dehning atoms or fragments in molecules. This 
has so far been investigated with Eq. (1) as the starting point. We refer the 
reader to Ref. [7] for a recent comprehensive review on density partitioning 
techniques, especially those aimed at dehning chemically and physically op¬ 
timal subsystems. Additionally, Refs. [11-13] contain outstanding reviews 
on the subject of subsystem DFT. 

In this topical review, we will touch on two competing concepts regard¬ 
ing the merits of subsystem DFT. They can be broadly summarized as the 
ability to: (1) reduce the computational scaling, and (2) gain a fundamen¬ 
tal understanding of the properties of the supersystem that otherwise would 
be hidden in the underlying complexity of the solution of a single electronic 
problem. For both of these points, we will present some of the most recent 
work (both published and unpublished) that our group has been producing 
in the past two years. These include: a new implementation of subsystem 
DFT in the plane wave basis set, an implementation of nuclear (ions) forces 
for Born-Oppenheimer dynamics, an implementation of real-time subsystem 
TD-DFT, and a formulation of a van der Waals theory based on subsystem 
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DFT. Points (1) and (2) are realized in each of the mentioned developments, 
although at different levels and with different outcomes. Each of the topics 
considered will be cast in the framework of the current state-of-the-art. 

This review is organized as follows. First, we will describe the two most 
common computational avenues to achieve a density partitioning, i.e. one 
employing kinetic energy functionals and the other one that does not. As our 
group focuses on the former, we will refer the reader to appropriate references 
for the latter when needed. We will also touch briefly on an alternative 
embedding trend, which partitions the wave function rather than the density. 
Second, we will scrutinize the underlying differences between the currently 
available implementations of subsystem DFT, such as basis sets and the 
choices of algorithms for solving the Self-Consistent Field (SCF) equations. 
After that, we will discuss the extension of subsystem DFT to the time 
domain by way of the Mosquera-Jensen-Wasserman theorem [14,15] (i.e., 
the extension to subsystem DFT of the Runge-Gross theorem [16]). This 
will provide us with a discussion of the linear-response and the real-time 
formalisms. 

1.2 Two flavors of subsystem DFT 

In this section, we will introduce two of the most popular methods to go 
about using Eq. (1). The hrst method makes use of the so-called nonadditive 
functionals to describe the interactions between subsystems. The second 
method, instead, requires the availability of the density of the supersystem 
and solves for those subsystem densities that reproduce the density of the 
supersystem exctly through Eq. (1). These two methods are very different in 
spirit, and here we will outline their main differences. 

It all starts from dehning the Lagrangian that needs to be minimized in 


7 



order to find the electron density of the supersystem, 


ll^DFT [p] = Ehk [p] + / next(r)p(r)dr - /i 


p(r)dr — iVe 


( 2 ) 


Herein, N^. is the number of electrons, fi the chemical potential, and Ehk [p] 
is the Hohenberg-Kohn functional [17] as reformulated by Levy [18] and 
Lieb [19]. The Kohn-Sham DFT (KS-DFT) reformulation of the above La- 
grangian [20] relies on the one-to-one mapping of Us-representable electron 
densities (we will only consider these type of electron densities in this review 
unless otherwise stated) on to a set of noninteracting electrons, also known as 
the Kohn-Sham system, whose wavefunction is described by a set of Kohn- 
Sham orbitals, {0*} (curly brackets indicate a set throughout this work). The 
orbitals are found by minimizing the following KS-DFT Lagrangian 

Eksdft [{0i}] = Ehk [p] + / next(r)p(r)dr - ^ eij [(0i|0j) - 5p]. (3) 

ij 

Usually, the above Lagrangian features a partitioning of the HK functional 
into 

Ehk [p] = Ts [p] -|- E^l [p] -|- E^c [p] , (4) 

where the electronic Coulomb repulsion energy (F^h), the noninteracting ki¬ 
netic energy (Tg) and exchange-correlation energy {E^c) functionals are in¬ 
troduced. 

By imposing Cksdft [{0i}] to be stationary, and only considering the 
so-called canonical solution (i.e., a diagonal eij matrix), the KS equations 



are recovered [5]. Namely, 


I 

--V^ + Wext(r) +WH[p](r) +t^xc[p](r) )0i(r) = €i(f)i{r). (5) 


Iks 


1.2.1 Use of nonadditive density functionals: The Frozen-Density 
Embedding method 

In the presence of multiple subsystems, each of the subsystem densities needs 
to integrate to a predetermined number of electrons, in order to ensure N- 
representability. For this purpose, we now consider the subsystem DFT La- 
grangian. 


[{P/}] — Ehk 

-Ns 

5^P/(r) 

+ / t'ext(r) 

-Ns 

X^P/(r) 

Ns 

j Pl(r)<ir - Nj 


. I 

J 

. I 

I 

(6) 


Cs [{P/}] is equivalent to Cdft [p] in Eq- (2) [5] with one redundancy — the 
total electron density is artihcially written as in Eq. (1) and the constraint of 
integration to the total number of electrons has been split into Ns constraints, 
one per subsystem. Of course, this partition has no effect on the result of 
the minimization, and having labeled the subsystem chemical potentials as 
fxi is hctitious at this stage, as the subsystem chemical potentials must all 
be equal to the supermolecular chemical potential [21-23]. 

A computationally practical avenue is found by rewriting Cs [{P/}] a 
slightly different way. Namely [24,25], 



-Ns 1 Ns 

5^P/(r) 


I 


I 


j p/(r)dr-iV7 

(7) 
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which once again is equivalent to Eq. (6). The nonadditive HK functional 
takes the form 

Eit [{p/}] = [{p,}]+[{p,}]+Er [{p/}] (8) 

where Eh, Tg, and E^c are not additive and yield the corresponding nonad¬ 
ditive functionals defined as, 

Ns 

F [p] = [p/, pn, • • •, Pats] + [pi ], with E = T^, E^c, and Er. (9) 

1 

Although the above formulation involving nonadditive functionals has 
been introduced only relatively recently [26-28], the idea of splitting a super¬ 
system into subsystems and approach the problem from a DFT prospective is 
much older and dates back to the early attempts in applying the “statistical 
model” [29-32] now more commonly referred to as orbital-free DFT [33,34], 

Minimization of the Lagrangian in Eq. (7) with respect to one subsystem 
density at a time (keeping the others frozen) yields a set of coupled differential 
equations. Namely, 

5Cfde [{pj}] _ 5Ehk[pi\ ^E^^ [{pj}] ^ 

(V,(r,.0) ipM (Jw(r).I)) (J„(r).0| 

( 10 ) 

In this review, the use of Eq. (10) is considered the definition of the Frozen- 
Density Embedding (FDE) method. Such semantical choice stresses the fact 
that the minimization of Cfde [{P/}] is carried out with partial functional 
derivatives as opposed to using full functional derivatives as proposed by 
Lahav and Kliiner [35]. In the density embedding literature, sometimes the 
FDE acronym is applied to a case where only one subsystem is included 
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in the variational problem Eq. (10) [28,36], and the term “subsystem DFT” 
is reserved for the self-consistent version (i.e. all subsystem densities enter 
the variational procedure). Throughout this review, we will use FDE and 
subsystem DFT interchangeably, as we view FDE [as dehned by Eq. (7-10)] 
as one possible algorithm to achieve a subsystem DFT treatment of the 
electronic problem. 

The dehnition in Eq. (10) can be cast in a KS-DFT scheme by transform¬ 
ing the Lagrangian in Eq. (7) to the following auxiliary KS-like Lagrangian 


Ns 


^KSCED [{.#.'}] = j2^.{p,]+Tr‘iiP,}]+ 

I 

Ns 

[pi] + [{P/}] + 

I 

Ns 

+J2E.C [pi] + [{pj}] + 


^'fxt(r)p/(r)dr + 


^eMpj{r) 

.J, K^I 


dr -I- 


Xl^ext(r)pj(r) 




dr 


Ni 


5^4 -%) 

Pi) I 


( 11 ) 


The external potential has been partitioned into subsystem contributions for 
convenience, i.e. Uext = Z)44xf In going from Cksdft [{0i}] to Cksced [{4}], 
the supersystem has been mapped onto a collection of auxiliary KS systems 
rather than a single one. Cksced [{4}] is labeled as Kohn-Sham with Con¬ 
strained Electron Density (KSCED) in line with the previous literature [37]. 

It is important to point out at this stage that we made no reference as 
to what dehnes the subsystems. All we have imposed is that the subsystems 
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must integrate to a preset number of electrons and be n^-representable. Thus, 
there is no unique solution to the variational problem involving Eq. (11). Dis¬ 
regarding for the moment the issue of non-uniqueness, imposing Cksdft [{0i}] 
to be stationary upon variation of a subsystem orbitals, {0f}, while keeping 
the other orbitals frozen, leads to the following one-electron equations which 
are given in two equivalent forms. 



where the nonadditive potentials are dehned as spifr)^^^ • 

We have also defined in Eq. (12) the embedding potential nemb(i‘) as a poten¬ 
tial collecting all terms linking the various subsystems. 

To summarize, in this section we have introduced the fundamental re¬ 
lationships at the base of the FDE method. In FDE, the term of central 
importance in the interaction between subsystems is the nonadditive kinetic 
energy functional (NAKE) whose exact expression is still unknown. Several 
approximations are available for it with a common starting point being non¬ 
interacting KE functional approximants, T, [33,38]. There is computational 
evidence [7,39,40] supporting the idea that the non-uniqueness issue is re¬ 
solved when employing approximate functionals as they effectively minimize 
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the inter-subsystem density overlap. 


1.2.2 Exact embedding: Partition DFT and Potential-Functional 
Embedding Theory 

Levy and Perdew [5,41] first, and then Wu and Yang [42,43] provided a mean 
to compute the KS effective potential when the input quantity is the electron 
density. This is achieved by first dehning the Wu-Yang functional 

t^WY [-^s] = Ts [p] + j ns(r) (p(r) - p(r)) dr, (14) 


which is to be considered a functional of the KS potential only, as the trial 
density, p, is directly obtained from the KS potential. Minimizing Cwy ["^s] 
w.r.t. the trial density p recovers the KS equations. Maximizing it w.r.t. 
the trial potential Vs (which acts as a continuous set of Lagrange multipliers) 
yields the unknown KS potential. Thus, the solution of the following coupled 
equations is sought; 


5Cwy ['^s] _ g 
5p{r) 

5Cwy [t^s] _ „ 
(5ns(r) 


-V^ -t-ns(r^ 


0i(r) 


ei0i(r). 


(15) 

(16) 


with the understanding that Eq. (15) involves a minimization, and Eq. (16) a 
maximization. A number of research groups, each with a different flavor of 
the theory, have followed the idea of reconstructing the embedding potential 
of FDE from the knowledge of the supersystem density [12,44-47]. The 
first numerical implementation of this technique is due to Roncero et al. [48] 
for regular subsystem DFT and Jacob et al. for a density partitioning with 
capping groups [49]. To aid our explanations, we subdivide the work of others 
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in this framework into two categories, both requiring the supersystem density 
as input. 

The first one prescribes a predetermined partition of the density so that 
Eq. (1) is achieved from the onset of the calculation [12,44,48-50]. Then, the 
unknown becomes the embedding potential that achieves such a partitioning. 
In the case of only two subsystems, this is sometimes cast in terms of a 
supersystem density, p, and a so-called “frozen” density, p//, automatically 
yielding p/ = p—pu provided that p(r) —p 77 (r) > 0 [22]. As the KE potential 
is equal to the negative of the KS effective potential plus a constant, the 
potential reconstruction technique is employed for recovering the effective 
KS potential for p/. 

The second flavor [7,9,46,47,51-53] poses the problem of finding the most 
appropriate partitioning of the system. In Partition DFT (PDFT), [7,9,47, 
51,52] this is done by dividing the total system into predetermined fragments, 
which can be either atoms, functional groups or separate molecules, and min¬ 
imize the sum of the fragment energies with the restraint that the sum of 
the fragment densities will match exactly the total molecular density. Given 
the constraint, the total molecular density is an input variable for the PDFT 
method and the total energy of the supersystem is never calculated. In po¬ 
tential functional embedding theory (PFET), the total energy is formulated 
using the embedding potential Vp. This leads to the following Lagrangian 


Ns 

^PFET [p, t'p] = [pi] + 


-Ns 

-p(r) 


dr. 


(17) 


which can be treated in the same way as Eq.(14), with Vp in PFET taking 
over the main role from Vg- It is quickly verified that Vp in Eq. (17) must 
be equivalent to Uemb of Eq. (12) provided that the embedding potential is 
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imposed to be the same for all subsystems [46]. We should mention that 
PFET and PDFT share Eq. (17) as the main working equation. In both 
theories the embedding (or partition) potential Up is unique and shared by 
all subsystems. 

The main drawback of these so-called “exact” embedding theories is that 
solving Eq. (15-16) is not computationally efficient. Despite this, the exact 
embedding method by Huang et al. was successfully applied to molecule 
surface interactions [11,54] in the context of wavefunction-in-DFT. Partic¬ 
ularly important is its application to the interaction energy surface of O 2 
with Al [54] - a long-standing problem in theoretical surface science. We will 
discuss the computational efficiency and effectiveness of the various flavors 
of subsystem DFT in section 2.1.1. 

1.3 Imposing Orthogonality: A revised Kohn—Sham 
scheme 

We have so far considered partitioning the supersystem by Eq. (1), e.g. tak¬ 
ing the electron density as the partitioned qnantity. An emerging trend is 
partitioning the snpersystem’s KS wavefnnction. That can be achieved by 
simply partitioning the KS orbitals into subsystems. Namely, 

{0j} -t |{0(i)7}, ■■ ■ , (18) 

The above partitioning leads to a trivial implementation of the method if 
all the orbitals are orthogonal to each other (and not only within a sub¬ 
system). The concept of imposing orthogonality originated in the held of 
psendopotentials, originally for the purpose of modelling core electrons and 
inclnding relativistic effects, but also for embedding [55-58]. An early exam- 
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pie of the benefits of imposing orthogonality between the orbitals of different 
subsystems comes from a post SCF calculation presented by Kolos et al. [59]. 
Upon orthogonalization, binding energies calculated with the “statistical 
model” [30-32] (otherwise characterized by low accuracy) were improved. 
Not imposing orbital orthogonality is inconvenient, as the Slater-Condon 
rules become somewhat more complicated [60,61]. In the KS method, or¬ 
thogonality is imposed in the Lagrangian, see Eq. (3). However, as mentioned 
before, one might want to extract information peculiar to only one portion of 
the supersystem and in order to achieve that from the onset of the calcula¬ 
tion, it is convenient to split the electronic problem into separate subsystem 
problems. 

In orbital space, such a partitioning can be achieved by the so-called 
Pauli blockade (PB) method [62] which was recently revived by Lukasz et 
al. [63,64], later by Miller and coworkers [65-67], and Gritsenko [22], as well 
as Hoffmann and coworkers [68,69]. The PB method relies on the effect of a 
penalty function to be added to the KS Fock operator introduced in Eq. (5), 
which reads 

Ns Nj 

Iks = Iks + X^7j X] 

Ji^i j=i 

where yj are small positive constants. 

It is important to remark that from Eq. (18) it is clear that the super¬ 
system is mapped onto a single KS system rather than to a collection of KS 
systems. This is an important point, as it prescribes omitting the NAKE 
potential in the working KS equations. Computational evidence of this can 
be found in Ref. [69] where the penalty in Eq. (19) was employed correctly 
(i.e. alone) and incorrectly (i.e. in conjuction with a NAKE functional). 

An overview of the different methods discussed in Sections 1.2.1, 1.2.2 
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and 1.3 are summarized in Table 1. 


2 Subsystem DFT for the ground state 

There is little doubt that solving the all-electron Schrodinger equation is 
computationally intractable. For example, if a general wavefunction is con¬ 
sidered (i.e., not necessarily a Slater determinant) the computational cost is 
(!1(A^!) with N being the size of the system (either the number of electrons 
or the number of basis functions employed). Avoiding the total solution of 
the Schrodinger equation while maintaining a satisfactory level of accuracy 
of the method is therefore the main goal of electronic structure theory. For 
example, if a Slater determinant alone is employed, the Hartree-Fock method 
is recovered with an associated formal computational complexity of (9(A^^). 
Unfortunately, this approximation has catastrophic consequences to the accu¬ 
racy of the method and nowadays no scientific publications use Hartree-Fock 
as the sole electronic structure method. 

DFT was formulated with the goal of hnding a method that would pro¬ 
vide electronic structures of molecular systems at a cheap computational 
cost while maintaining a satisfactory, or even exact level of accuracy. Since 
the Hohenberg and Kohn (HK) theorems were formulated [17], the electron 
density, p(r), has become the central quantity in electronic structure theory. 
This marked the birth of DFT. The HK theorems stated that once p(r) is 
known, then in principle also the molecular Hamiltonian of the system is 
available, and with that also the energy functional, E [pj. The vice versa is 
also achievable upon solving the appropriate Schrodinger equation. This is a 
vicious loop - as in reality neither the density nor the energy functional are 
known, and of course solving the full Schrodinger equation is computationally 
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intractable. 

The KS [20] self-consistent procedure has widened the scope of of DFT 
by hnding yet another direct relationship. This time between p(r) and an 
auxiliary system of noninteracting electrons whose wavefunction is a Slater 
determinant. Again, the Slater determinant is the gate to a computationally 
cheap method by avoiding tedious and expensive anti-symmetrization proce¬ 
dures. As for the HK theorems, KS theory gives no prescription regarding 
the energy functional which is still an unknown. In KS-DFT, the so-called 
exchange-correlation (XC) functional (i?xc [p]) needs to be approximated if 
practical calculations are to be carried out. Using local XC potentials yields 
a method that scales as 0{N^), and that generally is more accurate than HF. 
For this reasons, KS-DFT has become the standard method for determining 
electronic structures of molecule and materials. 

Because of the 0{N^) scaling with the system size, KS-DFT is limited 
in the system sizes it can approach. When modeling a physicochemical pro¬ 
cess, a careful choice of the model system precedes the computations. This 
is because, systems approachable by KS-DFT generally should not exceed 
500 atoms, so that the simulations can be completed in a reasonable time. 
However, it is now understood that if chemical accuracy in the predictions is 
sought, simulations must take into account the complexity of the environment 
surrounding the model system of interest [70-73]. The continuous struggle 
to increase the size of the systems in the full quantum chemical treatment 
is facilitated by the development of linear scaling techniques [74-77] and ef- 
hcient parallelization techniques [78,79]. In Section 2.1 we will show that 
subsystem DFT enjoys the beneht of being linear scaling (provided appro¬ 
priate approximations are made for the evaluation of the Hartree potential) 
and fully parallelizable in work and data. Section 2.2 will illustrate some 


19 



of the applications performed using FDE in recent years, demonstrating its 
additional advantage of offering physical insight into the nature of the super¬ 
system and the subsystems. 

2.1 Computational complexity and parallelization 

2.1.1 A divide-and-conquer approach? 

For the discussion of the computational speedup associated with subsystem 
DFT, we will focus on the implementations based on the methods discussed 
in Section 1.2.1, that involve the use of nonadditive density functionals. For 
this flavor of subsystem DFT, the computational gain can be achieved from 
two sources: (i) reduction of the dimension of the size of the problem that 
needs to be solved (i.e. one per subsystem) and (ii) the increased potential 
to parallelization. 

The original paper introducing FDE [28] proposed a simple concept: if 
the total density can be written as a sum of two densities p(r) = pi(r)-|-p 2 (i'), 
then it is sufficient to minimize the energy of the total system w.r.t. only one 
of the densities while keeping the other frozen to obtain the exact density 
and energy of the total system (provided that some mathematical conditions 
on the choices of the subsystem densities are satished). At a hrst glance, this 
appears to be an uncanny way to reduce the dimensionality of the problem 
as Nsub < Ntoti except for the following catch 22: one cannot impose the 
necessary conditions on the subsystem densities without the knowledge of the 
total density. As a result, one is obliged to perform the exercise iteratively, 
exchanging the roles of the frozen and variational density, also known as 
the freeze-and-thaw (FAT) cycles, until convergence is achieved. Of course, 
even then, one can argue that performing a calculation of dimension Ngub 
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several times is still more computationally advantageous than performing 
a calculation of dimension Ntot once because of the 0{N^) scaling, but for 
maximizing the computational potential of the method, additional techniques 
must be applied. 

The FAT procedure is designed to improve the initial guess of the frozen 
density, which needs to be such that at the end of the energy minimization, 
both subsystem densities (the frozen and the variationally minimized) are 
positive and n^-representable. Numerical tests have shown, however, that in 
certain cases the inclusion of the solvent as a frozen density without relax¬ 
ation through a FAT procedure, or only a partial relaxation withont reaching 
full convergence, might already provide satisfactory results [40]. For systems 
where the full relaxation of all snbsystems is necessary, FAT can be replaced 
by a mnch more computationally efficient procednre where all snbsystems 
are optimized simultaneously, with a total density update at each SCF itera¬ 
tion [80-82]. This is eqnivalent to constrncting a block-diagonal Fock-matrix, 
where each snbsystem is represented by a block. Both approaches lead to 
the same minimum if brought to a total point of convergence, where the for¬ 
mer option might be performed partially and the latter option holds a much 
stronger potential for parallelization, as will be discussed below. 

When working with the linear combination of atomic orbitals (LCAO) 
method for the molecular orbitals, the 0{N^) scaling is determined by the 
number of basis set functions used in the calculation. A straightforward 
implementation of FDE would involve the use of the supersystem basis set, 
which would not resnlt in computational speed np bnt can be usefnl for per¬ 
forming benchmarking calculations [83,84]. A more compntationally efficient 
option is to use a basis set centered only on the atoms of the snbsystem in 
question. Although, depending on the nature of the system being stndied 
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it might be necessary to include the orbitals of other subsystems as well 
in order to have a basis set flexible enough to model the flow of electron 
density from one fragment to another [85]. Additionally, one can opt for 
calculating the density of the frozen system using a htting set instead of 
using the exact density [86] and employing efficient numerical integration 
schemes [81,87] for large environments. In the subsystem DFT implementa¬ 
tion in ADF [86], the numerical integration grid is centered on the nonfrozen 
subsystem and therefore, for sufficiently large environments, the integration 
grid becomes independent of the size of the environment and, as a result, 
so does the calculation time of the numerical integration. Another example 
where subsystem DFT is combined with a numerical integration scheme using 
hierarchical real-space grids is by Shimojo et ah [81], where, combined with a 
high level of parallelization, it achieves linear scaling and a parallel efficiency 
of 0.985 on 128 processors. The flexible implementation of subsystem DFT 
in ADF [86] allows to choose which subsystems will be fully nonfrozen, which 
will be frozen but relaxed using FAT cycles to full or partial convergence, and 
which will be kept completely frozen (Fig. 1). Such a set up is, for example, 
of use when studying a molecule in a solvent, where the molecule will be the 
nonfrozen fragment, the hrst solvation shell relaxed frozen and the rest of 
the solvent molecules fully frozen. 

[Figure 1 about here.] 

Alternatively, subsystem DFT can be implemented in a code for periodic 
systems using plane wave (PW) basis sets [82] which are not centered on 
atoms. Such codes have very distinct advantages for applications on periodic 
systems and condensed phases and have a high potential for parallelization 
due to the use of a 3D space grid (See for example [78,81,88]). Reducing 
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the dimensionality of the problem in subsystem DFT calculations in PW- 
based codes is, however, more challenging since the basis set size depends 
on the size of the periodic cell, which, in a straightforward implementation, 
is determined by the size of the supersystem. Reducing the dimensionality 
in this case is more involved than in codes using localized basis sets, but 
not impossible. One can implement the use of smaller interlocking cells for 
each subsystem, cut out of a larger supersystem cell. In this case, the KS 
orbitals of each subsystem can be expanded using the plane waves spanning 
the smaller cells, while the density and the potential are evaluated using the 
plane waves spanning the supersystem cell. This concept is illustrated in Fig. 
2 . 


[Figure 2 about here. 


2.1.2 Subsystem DFT: a naturally parallel method 

The efficiency of the subsystem DFT implementation is strongly influenced 
by the efficiency of the quantum chemical code in which it is implemented. 
Any advanced computational techniques such as linear scaling, numerical 
integration, nearsightedness and parallelization that are present for the so¬ 
lution of the supersystem KS-DFT problem will be transferable to the sub¬ 
system DFT problem. Subsystem DFT offers above these optimizations a 
higher level of parallelization for both work and data. This concept is illus¬ 
trated in Fig. 3a using our implementation of subsystem DFT in Quantum 
ESPRESSO [82,89]: Each subsystem calculation can be seen as a separate 
KS-DFT calculation optimizing the KS orbitals of the subsystem. In this 
implementation, the standard option is to update the density between the 
subsystems after each SCF cycle, though a FAT option also exists. Needless 
to say, both possibilities lead to the same numerical results. At the beginning 
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of each SCF cycle, the total electron density (pin(i‘)), composed of the sum 
of the subsystem densities, is communicated to all the subsystems for the 
evaluation of the potential. At the end of each SCF cycle, the new subsys¬ 
tem densities are reduced to form a new total electron density (pout(i')). 
Such scheme is easily implemented using MPI, but one must take care to 
balance the computation time required by different subsystems to perform a 
single SCF iteration. In the Quantum ESPRESSO code this is achieved 
by optimizing the number of processors assigned to each subsystem when 
the subsystems vary in size. During the SCF cycle, the processors assigned 
to each subsystem are grouped together in an “intradmage” communicator 
and the code uses the standard parallelization of Quantum ESPRESSO, 
where the 3D real and reciprocal space grids are distributed between differ¬ 
ent processors. Between the SCF cycles, the density of each subsystem is 
reduced to an array only allocated on the head process of the communicator, 
and subsequently reduced and distributed among the subsystems using the 
“interJragment” communicator. This concept is illustrated in Fig. 3b. 

[Figure 3 about here.] 

2.1.3 Non-additive functional approximants 

We have seen in Section 1.2.1 that one way to cast DFT in a subsystem 
fashion involves the use of nonadditive functionals, i.e., each energy term of 
the supersystem can be expressed as the sum of additive and nonadditive 
contributions. From Eq. (9) we have 

Ns 

F[p] = 5^F[p/] + [{pQj, with F = Ts, E,,, and Eh (20) 
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where the nonadditive term is 


Ng 

f’“‘‘lW]=f W-E^|P/]. (21) 

I 

While the above equation is almost trivial for the Hartree and semilocal 
XC terms, problems arise in the evaluation of the noninteracting KE of the 
supersystem, when the total density is built from the KS orbitals of the 
subsystems, which are not required to be orthogonal to each other across the 
fragments. 

It is common practice to evaluate the NAKE using pure functionals of the 
electron density rather than the (unknown) KS orbitals of the supersystem; 

Ng 

r;™" [{p/}] = % Ip] - Y7, Ip,I (22) 

I 

where T), is an approximate functional for the noninteracting KE. 

Clearly, the accuracy of the FDE results will be strictly related to the 
quality of the approximate NAKE functional employed in the calculation. 
Over the past decades many KE functionals have been developed. As in 
the case of the exchange-correlation functionals, there are different types of 
KE functional approximants, spanning from the local TF/LDA functional 
[90,91], to the semi-local GGA family [92-99], to fully nonlocal functionals 
[33,34,100]. 

Since the only approximation introduced in the formulation of FDE is the 
use of approximate kinetic energy functionals for the evaluation of the non¬ 
additive energy (and potential), it follows that a systematic way to expand 
the types of systems this method can effectively simulate is to improve the 
KE functionals. We refer the reader to Ref. [12] for a comprehensive review 
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on (among others) the NAKE functionals. 


2.2 Applications of subsystem DFT 

2.2.1 Molecular systems: localized basis sets 

Since its birth [28], the FDE method has been implemented in a wealth of 
Quantum chemical packages. Examples include Amsterdam Density-Functional 
(ADF) [101], TURBOMOLE [99,102], deMon [36,103], Dirac [104], Q-Chem 
[45], and MolCas [105]. The number of features and the flexibility of the 
implementations may vary, but all of them share some important similari¬ 
ties, such as the use of localized atom-centered basis sets (Slater or Gaussian 
type). Where available, the FAT scheme [36] to achieve self-consistency is 
employed. 

The FDE method combined with the currently available semilocal NAKE 
functionals has been proven very successful in predicting the electronic ground- 
state properties of a vast array of noncovalently bonded multi-molecular 
systems [101,106-111]. Examples include an accurate prediction of ligand 
dissociation and proton transfer in metal complexes [85,112], the correct 
interaction energy and electronic structure for weak and strong hydrogen¬ 
bonding systems [113,114], van der Waals interactions [84,115], and up to a 
certain extent also reproduced the right electronic structure of Lewis acid- 
base complexes characterized by weak dative bonds [116]. 

[Figure 4 about here.] 

A powerful flavor of subsystem DFT is 3-FDE. It has been proposed by 
Jacob and Visscher [49,117] to simulate peptide chains and potentially full 
proteins. Starting from the Molecular Fractionation with Conjugate Caps 
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(MFCC) scheme [76,118,119], the chain is cut at the peptide bond, and 
capping groups are added at the extremes to mimic the original environment. 
The total density of the system is then the sum of the capped amino acids 
minus the density of the capping groups (See Fig. 4). 3-FDE introduced 
several novelties. The densities are not just those from isolated fragments (as 
in MFCC), but are calculated using a FDE embedding potential, and possibly 
converged using freeze-and-thaw, it can be applied to subsystems connected 
by covalent bonds, whereas regular FDE would fail due to the approximate 
NAKE functionals used. 3-FDE has been applied to several peptide pairs and 
to the ubuquitin protein, yielding signihcantly more accurate results than 
the MECC scheme. More recently, it has been extended to the calculation 
of excitation energies. [117,120] 

2.2.2 Condensed-phase systems: plane wave basis sets 

Subsystem DFT implementations for periodic systems do not have a history 
as extensive as the molecular codes counterpart, mostly because the cur¬ 
rent existing implementations (such as CP2K [80,121], ABINIT [122,123], 
CASTEP [35]) have limitations of different nature, by either limiting the 
number of fragments to 2, or by not having the possibility to sample the hrst 
Brillouin zone (FBZ). 

Efforts in our group to £11 this gap are ongoing [82] and involve the 
implementation of FDE in the plane wave (PW) code Quantum-Espresso 
(QE) [89,124]. Using PW allows for a very accurate and efficient calculation 
of the Hartree energies and potentials in reciprocal space. Also the evaluation 
of the forces acting on the nuclei is more straightforward compared to a 
localized basis set implementation (see Sections 3 for details and application 
to molecular dynamics). Moreover, using a PW basis makes the system 
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intrinsically periodic, making it the most natural basis set to expand the 
Bloch states of such systems, thus allowing for a sampling of the FBZ with 
arbitrary accuracy. Finally, in the subsystem DFT case, we have the option 
of simulating the fragments using different accuracies in the sampling of the 
FBZ, potentially leading to additional computational time savings. 

A disadvantage of a PW implementation is that the use of pseudopoten¬ 
tials in place of the nuclei and core electrons is always needed in order to 
avoid expanding the one-electron wavefunctions in the fast oscillating nodal 
region close to the nuclei. Specihc to subsystem DFT is also the fact that if 
the same simulation cell is employed for all the fragments, we end up having 
to solve Ns diagonalization problems in the same (large) PW basis set, while 
localized orbitals implementations have the ability of only including the basis 
functions centered on specihc fragments. This means that a speed up in the 
calculation through a divide and conquer approach does not come as easy 
for a PW implementation as it does for a molecular code (See Section 2.1.1). 

QE offers an ideal platform for the implementation of subsystem DFT, 
which enables the application of the method to solid state fragments such as 
surfaces. [82,89] The FDE implementation is added as a higher level of par- 
allellization to the existing code using the popular Message Passing Interface 
(MPI) libraries. To the best of our knowledge, ours is the most compre¬ 
hensive implementation of subsystem DFT in a PW basis. Specihcally, the 
implementation features simultaneously: (i) similarly to CP2K and most 
molecular codes, such as ADF, the ability of simulating an arbitrary number 
of fragments yielding self consistent electron densities, (ii) a fragment specihc 
sampling of the FBZ (CP2K only samples the F point, while the other peri¬ 
odic implementations are restricted to only model two subsystems at a time), 
and (iii) similarly to CP2K, calculation of the energy gradients wrt nuclear 



displacements. To achieve full self consistency, our implementation does not 
use the FAT procedure as a default. Instead, thanks to our MPI machin¬ 
ery and similarly to the CP2K implementation, we are able to overcome the 
active-frozen fragment distinction, simulating all the subsystems at the same 
time. The SCF iterations are synchronized across the fragments and a new 
embedding potential for each fragment is calculated at each iterative step. 
In addition, our implementation allows us to allocate an arbitrary number of 
processors to each subsystem depending on their size (see Section 2.1.2). 

The formulation of the theory for a periodic implementation of subsystem 
DFT is very similar to that of the molecular case as it has been explained 
in Section 1.2.1. However, important distinctions include the fact that we 
might have fractional occupations of the KS orbitals (either to help SCF 
convergence or to simulate an electron hnite temperature), and that the 
sampling of the FEZ has to be taken into account in the equations. As a 
consequence, the Janak’s kinetic energy, Tj, of a fragment is found in place 
of the noninteracting kinetic energy, T^, 
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where the subsystem Bloch waves are 0f(r) = e**^''ufj^(r). The electron den¬ 
sity of a fragment with respect to the KS orbitals is given by 

P/(r) = ^ / dk iMikWT- (24) 

“BZ J-BZ , 


[Figure 5 about here.] 

We proceeded [82] to benchmark our code against systems that had al¬ 
ready been studied using localized orbitals codes, observing very similar re- 


29 



suits: a good agreement with the KS reference for H-bonded systems, and a 
not completely satisfying agreement for Lewis acid-base complexes. We then 
proceeded to study for the hrst time the interaction of molecules and surfaces 
using subsystem DFT. Test systems included methane on a Pt(lOO) surface 
in eight possible conhgurations, and a liquid water bilayer on a Pt(lll) sur¬ 
face (a system comprised of 13 fragments). We showed [82] that subsystem 
DFT carried out with GGA NAKE functionals yields good results when the 
interaction between the molecule and the surface is dominated by electro¬ 
statics, like for instance the system shown if Figure 5c. On the other hand, 
due to the known shortcomings of the available GGA NAKE functionals, 
the quality of the results quickly degrades as the nature of the interaction 
becomes covalent, as it was observed in only two of the possible methane on 
Pt conhgurations (Fig. 5a) and by a lesser extent for a parallel water on a 
Pt(lll) surface (Fig. 5b). 

An important conclusion of our preliminary work is that in employing 
NAKE functionals with increasing accuracy, the FDE modeling becomes 
more accurate. 

2.2.3 Note on self-interaction 

It has been shown [125-127] that employing hybrid XG functionals (i.e. func¬ 
tionals including a fraction of the HF exact exchange) for the evaluation of the 
intra-fragment XG energy, greatly improves the quality of FDE for charge- 
transfer complexes, w.r.t. both the interaction energy and the predicted self 
consistent electron density. This can be rationalized by the fact that hybrid 
functionals are able to reduce the DFT self-interaction error arising from the 
use of approximate exchange functionals, effectively localizing the electron 
densities of the fragments. When the electron densities of the fragments 
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are less diffused, also the overlap between densities of different fragments is 
reduced. 

The subsystem DFT method has been exploited to constrain spin and 
excess electron densities on specihc fragments [127-131] even when KS-DFT 
fails to do so due to self-interaction error. This behavior of FDE arises 
from the fact that the FDE theory does not impose orthogonality between 
the orbitals of different subsystems, thus the typical delocalization of the KS 
orbitals originating from the orbital’s hybridization is completely avoided. In 
addition, and probably more importantly, spurious repulsive potential walls 
arise between the fragments due to the employed semilocal NAKE functionals 
at the location of the atomic shells of the frozen subsystems [44,128,131]. 
It is, therefore, possible to use FDE to approximate the wave function of 
diabatic states in processes involving charge transfer between a donor and 
an acceptor fragment [126,128,129], potentially mediated by bridge molecules 
which are also treated in a subsystem DFT fashion [130]. This framework 
has recently been applied [130,131] with great success to systems of biological 
interest and realistic size (such as portions of the DNA double helix), yielding 
qualitative and quantitative agreement with experimental data. 

It is important to point out that the diabatic states generated with FDE 
can only feature charge localization on an entire subsystem. If the diabatic 
state sought is to feature a charge localization on only a portion of the sub¬ 
system, then techniques such as constrained DFT [132] should be coupled 
with FDE. 


[Figure 6 about here.] 

With our plane wave FDE implementation, we have also ventured into the 
possibility of localizing charge and spin onto a single subsystem. There 
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are some processes in surface science and condensed phase that are not 
straightforwardly modeled with KS-DFT due to the self-interaction error 
introduced by GGA XG approximants. Two cases are particularly interest¬ 
ing, small molecules interacting with metal surfaces, and solvated radical 
species. When metal surfaces are involved, the simplest reactive scattering 
experiments, such as the ones involving diatomics, have in recent years be¬ 
come theoretical chemistry puzzles [most notably 02-l-Al(lll)] because of the 
complexity introduced by the self-interaction [54,133-136]. With the novel 
FDE implementation, our group has performed calculations on HO* embed¬ 
ded in a water solvent. In Figure 6 a snapshot of the spin-density extracted 
from an ab-initio molecular dynamics simulation of OH radical embedded in 
11 closed-shell water molecules shows that the spin-constrained FDE spin- 
density is correctly localized on the OH fragment reproducing experimental 
observations [137] and self-interaction corrected KS-DFT simulations [138]. 
When a semilocal KS-DFT calculation is carried out, the spin-density is too 
delocalized [138,139] and does not compare favorably with the experiment. 

3 Born—Oppenheimer molecular dynamics 

A very popular approach when tackling molecular dynamics (MD) is to rep¬ 
resent the coupled electron-nuclei time-dependent wave function through the 
adiabatic approximation 

OO 

{R/}; t) = ^n({K}; {R/} Xn({R/};t) ^ Tfc({ri}; {R/} Xfc({R/};t) 

n=0 

(25) 

where represents the n-th solution of the electronic Hamiltonian in the 
clamped-nuclei framework, and Xn are time dependent nuclear wavefunctions. 
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Representing the total wave function as a product of uncoupled nuclear 
and electronic wave functions allows us to describe the motion of the nuclei 
by solving a time-dependent Schrodinger equation involving only the nuclear 
coordinates, known as the Born-Oppenheimer approximation: 

where is the potential energy surface (PES) for a given k-ih elec¬ 

tronic state. Directly solving such equation is not a viable option for most 
systems, as among other things it would require the knowledge of a (3iV — 
6) dimensional PES. A popular formulation of ab-initio MD is the Born- 
Oppenheimer MD (BOMD), where the trajectory of the nuclei is propagated 
classically on a PES calculated on the fly from a ground state electronic 
structure calculation 


.0 
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The forces acting on a nucleus are then calculated using the Hellmann- 
Feynman theorem, and including Pulay force corrections due to the incom¬ 
pleteness basis set used (in the case of an atom centered basis). 

If a basis set independent of the position of the nuclei is employed (as 
plane waves are for instance), there are no Pulay forces. However, in most 
plane waves calculations pseudopotentials (PP) are used to mimic the nuclei 
and core electrons potential. Some pseudopotentials (e.g. ultrasoft [140-142] 
and PAW [143]) contain terms that are determined self-consistently and, at 
the same time, depend on the positions of the atoms. In those cases, Pulay 
corrections to the forces need to be included. Most modern PP are split 
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in a local radial part (needed to represent electrostatics) and in a nonlocal 
term (i.e. a term depending on the angular momentum of the single particle 
orbital it is applied to). In a subsystem DFT calculation employing PP, the 
nonlocal part of the PP of atoms belonging to one fragment only apply to 
the orbitals of that particular fragment, while the local part is shared across 
the subsystems [80,82], This ensures that the core electrons are assigned to 
only one auxiliary KS subsystem. 

There are not many examples of BOMD based on a subsystem DFT elec¬ 
tronic calculations. The CP2K implementation of FDE was recently used to 
model dynamical molecular properties of solvated molecules, with particular 
focus on the dipole moment [144]. One study from Kodak and Bernholc [145] 
proposed a method combining KS-DFT and a frozen density environment to 
simulate solvated biological systems. In that study, the periodic simulation 
cell is split in two spatial domains, see Figure 7. Atoms belonging to a 
smaller inner cell (typically a molecule and its first solvation shell) are simu¬ 
lated using a full KS method, while the solvent molecules outside this region 
are modeled by a hxed parametrized electron density: A linear combination 
of Gaussian functions, with one Gaussian function centered on each atom. 
Along the dynamics simulation atoms can move in and out of the fully quan¬ 
tum mechanical cell. Despite the fact that this method is not strictly speak¬ 
ing subsystem DFT, the embedding potential of the environment on the QM 
region is the same as Eq.(13) (e.g. it includes the nonadditive potentials). 

[Eigure 7 about here.] 

Molecular dynamics of liquid water using FDE have also been studied 
by lannuzzi et al. [80]. FDE had already been proved capable of predicting 
the optimal geometry for the water dimer [107], but had not been tested 
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on liquid water. lannuzzi et al. [80] carried out BOMD simulations for 64 
independent water molecules, using an array of NAKE functionals. It was 
shown that none of the NAKE functionals available at the time allowed FDE 
to correctly reproduce the experimental radial pair distribution functions, see 
Figure 8. Most notably FDE would fail to reproduce even qualitatively the 
second solvation shell in the 0-0 distribution, go-oi'f')- A posteriori, this 
result might come as a surprise, because in 2012 Hu et al. [146] showed that 
by using an approximate effective embedding potential htted from dimer 
data, a qualitative agreement with the experimental radial distribution of 
liquid water was achieved. 

[Figure 8 about here.] 

The ability of running molecular dynamics and geometry optimization 
has also recently been included in our periodic subsystem DFT implementa¬ 
tion [82,89]. Coupled with the ability of the code to sample the FEZ of the 
fragments, the code provides a platform for simulating for example the dy¬ 
namics of molecules on surface taking advantage of subsystem DFT. Current 
efforts in our group are directed towards applying the code on the dynamics 
of condensed systems, including shedding light on the elusive liquid water 
simulations. 


4 Time-dependent subsystem DFT 

4.1 Extension of the Runge—Gross and van Leeuwen 
theorems 

In the ground state, the Hohenberg-Kohn theorem [17] maps uniquely the 
wave function to a n-representable density, by proving that two potentials 
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differing from each other by more than a constant cannot produce the same 
density. For the time-dependent case, a similar mapping is done by the 
Runge- Gross theorem. [16] The full quantum chemical description of a time- 
dependent problem is given by the time-dependent Schrodinger equation 

( 28 ) 

where the time-dependent Hamiltonian is partitioned into a static part and 
a time-dependent part, which is a direct consequence of the application of 
an external time-dependent potential, and could be written as: 

H = Ho + Hi{t) (29) 

This connects a certain time-dependent external potential n(r, t) to the time- 
dependent wave function ^(t), given an initial state Tq. As in the ground 
state, a potential v'{r,t) that differs from v{r,t) by more than an additive 
time-dependent scalar function c{t) will result in a different wave function. 
The Runge-Gross theorem [16,147] proves that there is an analogical cor¬ 
respondence between the time-dependent potential and the time-dependent 
density p{r,t). The situation is however slightly more complicated than in 
the case of the ground state due to the presence of an initial boundary con¬ 
dition, namely the initial state Tq. In order to prove the one-to-one mapping 
between the time-dependent potential and the time-dependent density, one 
needs to resort to the concept of current density Without going into 

the mathematical details, one can show that, when starting from a certain 
state and two different potentials, the two current densities will start diverg¬ 
ing from each other at a time inhnitesimally larger than to and, as a result, 
also the corresponding time-dependent densities. In other words, the time- 
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dependent density is fully determined by the time-dependent potential and 
the initial state. This also implies that the time-dependent Hamiltonian and 
wave function are functionals of the time-dependent density. 

The van Leeuwen theorem [148-150] can be seen as the Kohn-Sham the¬ 
orem equivalent of TDDFT, since it allows us to link a fully interacting sys¬ 
tem to a noninteracting system, resulting in the well known time-dependent 
Kohn-Sham (TDKS) equations 

[-^V^ + n,(r,t)]0i(r,t) = (30) 

where instead of evolving the full density in time, one evolves the one-particle 
Kohn-Sham orbitals. We know from the Runge-Gross theorem that a certain 
time-dependent density p(r,t), given a certain initial state Tq, corresponds 
to an external potential v{r,t). This system will have a particular electron- 
electron interaction t(;(|r — r'|). van Leeuwen’s theorem states that there 
will exist a different system featuring a different electron-electron interac¬ 
tion tc'dr — r'l), associated with a different external potential v'{r,t) and a 
different initial state <ho that will be associated with the same time depen¬ 
dent density. In the special case of w'{r,t) = 0, this is the noninteracting 
Kohn-Sham system. 

Both the Runge-Gross and the van Leeuwen theorems have several prac¬ 
tical restriction [151] such as the Taylor expandability of the potential and 
the density, as well as the n-representability of the density. The details of 
these conditions and their breaking point he outside the scope of this re¬ 
view. The relevant question here is how the two theorems can be viewed 
in the framework of the different flavors of subsystem DFT. As discussed in 
Section 1.2.2, there are two flavors of subsystem DFT, centered around non- 
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additive density functionals and around unique embedding potentials. The 
formulation of subsystem DFT influences the way the Runge-Gross and van 
Leeuwen theorems should be viewed. In Frozen Density Embedding (FDE), 
the energy is minimized w.r.t. the density of each subsystem, each of them 
representing a separate KS system, and therefore also in the TDDFT ex¬ 
tension, the Runge-Gross and van Leeuwen theorems have to be discussed 
on a subsystem level. On the other hand, potential functional embedding 
theory (PFET) and partition density functional theory (PDFT) minimize 
the energy not w.r.t. the subsystem density, but with respect to an embed¬ 
ding potential, which is constrained to be a global quantity that is shared 
by all subsystems. Therefore, the Runge-Gross theorem for these methods is 
formulated with respect to the global embedding potential. 

At the beginning of a time-dependent FDE calculation, [152-154] each 
subsystem represents a Kohn-Sham system, mapped, on the one hand, to a 
noninteracting single Slater determinant wave function and an effective 
potential nf(r), and on the other hand, a many body wave function and 
an external potential v\r). As discussed by Gritsenko [22], when using 
approximate kinetic energy functionals for the nonadditive kinetic energy, the 
external (or effective) potential of each subsystem within the FDE method 
differs from the external (or effective) potential of the the supersystem by 
the error made by Ts[p]. For each subsystem I, the external (or effective) 
potential becomes 
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where ATs[p] refers to the error made by the approximate kinetic energy 
functional Ts[p] with respect to the exact kinetic energy Ts[p]. For the pur- 
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pose of the Runge-Gross theorem, we will first concentrate on the tl'^ as the 
initial state of each subsystem and the time-dependent external potential 
The question that rises is whether one can perform the proof of 
the Runge-Gross theorem using the subsystem external potential to show a 
unique mapping between t) and p^(r, t). The proof of the Runge Gross 

theorem relies on one major condition, namely the Taylor expansion of the 
potential about the initial time. Since the subsystem potential differs from 
the supermolecular potential by the kinetic energy functional error and does 
not exhibit singularities in time, one can assume this condition holds for all 
subsystems which are part of a supersystem, in which the potential is also 
Taylor series expandable. For the van Leeuwen theorem, an additional con¬ 
dition is necessary: the density must be analytic in time at to- Usually, it 
is assumed that when starting from a system that is in the ground state at 
t < to; this condition is satisfied. However, as shown by Maitra et ah [155], 
densities can nonetheless become nonanalytic in time in several special sit¬ 
uations. Since all subsystems are in the ground state at t < to, we for now 
assume that this condition is satished, but further research regarding the 
influence of the density partitioning on the analytic time behavior of the 
subsystem densities is needed. 

An important result is the unique mapping between the subsystem time- 
dependent potential t) and the subsystem time-dependent density p^(r, t). 

Since the time-dependent potential can be expressed as the sum of the iso¬ 
lated subsystem potential plus and embedding potential [Eq. (12)], there is 
also a unique mapping between the time-dependent density and the time- 
dependent embedding potential. In other words, the embedding dehnes the 
time evolution of the subsystem in a unique way. 

The formulation of the Runge-Gross theorem in the case of the time- 
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dependent PFET [15] and PDFT [14] theories is different, since there the time 
evolution is determined by the initial state 'Fq and the embedding potential 
'yemb(i', t). In ground state PFET, there is a unique mapping between the total 
electron density and the embedding potential, as opposed to the external 
potential in DFT. The Runge-Gross theorem has been recast using the time- 
dependent embedding potential hrst by Mosequra et ah [14] in the framework 
of the fragment-based PDFT [14] method, and later by Huang et ah [15] in 
the framework of the TD-PFET method. Both proofs show that if a system, 
dehned at t = 0 by an initial state and an embedding potential nemb(i', ^ = 0), 
evolves into two different time-dependent potentials nemb(i‘G) and 
that yield the same time-dependent total electron density, the two embedding 
potentials can only differ by a time-dependent scalar function. The Runge- 
Gross theorem can therefore be reformulated using the embedding potential. 

In practical calculations, the time evolution can be solved either through a 
linear response mechanism, which will be discussed in Section 4.2, or through 
a direct integration in time of the time dependent Kohn-Sham equations [Eq. 
(30)], which will be discussed in Section 4.3. 

4.2 Linear response approximation 

The ground state and excited states of a system which do not evolve with time 
are solutions to the time-independent Schrodinger equation. In principle, any 
ground state formalism can be extended easily to obtain excitation energies 
by considering slight deviations around the ground state, assuming that the 
correlation effects dominant in the ground state will also be mostly dominant 
in the excited state. The search for such a formalism led to the development 
of the Linear Response (LR) method that describes the response of to an 
applied weak perturbation of a system in its stationary state. Since the 
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perturbation is weak, the response can be truncated at the linear order to a 
good approximation. The LR method can be applied as an extension to any 
ground state method for energy calculation and when extended in the context 
of DFT in a time-dependent regime, is known as the Linear Response Time- 
Dependent DFT (LR-TDDFT). Since its formulation, it has grown to be 
very popular and have been applied to chemical problems of varying degrees 
of complexity [147,156,157]. 

LR-TDDFT allows to obtain excitation energies by performing the fol¬ 
lowing gedanken experiment. Starting from a static system in an initial 
state, one applies a time-dependent external potential of the form Uea;t(i', t') = 
—t) and monitors the response of the system. Uo(i') is the 
nuclear potential to which the system is subjected before t and Vappi{r,t') is 
the explicitly time-dependent perturbation turned on at time t using the step 
function 6. This external potential induces a nonstationary state 

which is described by a linear combination of ground and excited states. 
Hence it should be possible to obtain the excitation energies of the system 
from the structures of the nonstationary wavefunction. 

As discussed in Section 4.1, for each time-dependent density, there is 
one to one mapping between the fully interacting system evolving under the 
influence of Uext(i',^) and a KS system evolving under the influence of an 
effective potential Veff{r,t). The external (or effective) potential at different 
times t' and positions r' couples with the density of the system, causing 
density changes within the system. For a weak external perturbation, the 
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density change is linear in the potential and the linear response is given by: 


(5p(r,t) 



dr'xir, r', t - t')5vappi{Y', t') 
dr'xo(r, r', t - t') 


(32) 


where x(r,r',t — t') is the response fnnction of the fnlly interacting system 
and ~ t') is the response fnnction of the KS system. Now, staring 

from the von Nenmann eqnation, given by: 

^^^P=[H,P] (33) 


in the Heisenberg representation [158] one obtains the response fnnction in 
the time domain: 


X{r,r',t-t') = -ie{t-t'){'^o\[p{r,t),p{r',t')]\'^o) (34) 

where 6 is the step fnnction. Applying a Fonrier transformation into the fre- 
qnency and inserting the resolntion-of-identity, the Lehmann representation 
of the freqnency dependent density-density response [151] is obtained: 

y f (^o|p(i-)|^n)(4/.|p(r0|4/o) (d/o|p(r0|4/.)(vl/,.|p(r)|vl/o) -| 

’ ’ ^ I U — iln + ip U + iln + ip J 

n=l ' ' 

(35) 

where It^o) is the gronnd state wave-fnnction of the many-body system and 
= En — Eo is the excitation energy. Eq. (35) allows us to extract excitation 
energies directly from its pole structure: it will have poles whenever the 
frequency of the probing potential matches with the corresponding excitation 
(or deexcitation) energy of a particular excited state n. However, owing to the 
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time lag in the application of the external perturbation and the response of 
the system, a small positive shift rj appears in the denominator in the complex 
region of the frequency plane [151]. Thus, the poles of the equation he in the 
upper region of the complex frequency plane and its strength depends on the 
value of the numerator (transition densities). 

As x(r, r',oj) is the density-density response function of a fully interacting 
system, the poles of Eq. (35) have the correct structure corresponding to 
all states. The KS equivalent of Eq.(35) involves the use of single Slater 
determinants, which leads to: 


Xo(r,r',ca) = 

3,k 


0fc(r)0i(r)0;(rO0fc(E) 

u - {Sj - Ek) + irj 


(36) 


where fk and fj are occupation numbers of orbitals k and j, respectively. 
For xo(r,r',a;), the non-zero contributions are only when j is an occupied 
orbital and k is an unoccupied orbital, since all other terms cancel out. 
The KS response function only has poles whenever the energy difference 
{ej — Ek) matches the frequency of the applied external perturbation. Com¬ 
paring x(r,r',a;) and xo(r, r',a;) reveals an important difference: The true 
density-density response function in Eq. (35) has multiple solutions of all pos¬ 
sible excitation rank while for the KS density-density response function, the 
solutions in Eq. (35) are limited to single excitations. As a result, Xo(i', r',uj) 
cannot reproduce the correct pole structure of the fully interacting system. 
The connection between xo(I■^J'^|^) and x(r,r',a;) is readily extracted from 
Eq. (32): 


X ^(r,r',a;) = Xo^(r,r',a;) 



/^c(r,r',a;) 


(37) 
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where f^c is the frequency-dependent exchange-correlation kernel , given by: 




6vxc[p\{r,t) 

6p{r',t') 


[po]' 


(38) 


The frequency-dependent exchange correlation kernel is the key variable 
which would generate additional solutions beyond the KS response function 
and give correct structnres and nnmber of the poles for a true interacting sys¬ 
tem. However, since the exact exchange correlation fnnctional is nnknown, 
approximations must be introduced for practical applications. Among them, 
there is the adiabatic approximation which involves neglecting the freqnency 
dependence of the exchange correlation kernel, or in other words, employing 
a static exchange correlation potential which has no dependence on the time- 
evolntion of the entire system at different time t' and position r'. The effect 
of eliminating the freqnency dependence is readily seen from Eq. (36): any 
freqnency independent shift in the denominator cannot generate additional 
solntions bnt can only shift the positions of the poles. Hence, the possibility 
of acconnting of excited states rich in mnltiple excitations is forgone in the 
adiabatic approximation. Additionally, the exchange correlation kernel be¬ 
comes local in natnre with repercussions to the description of states nonlocal 
in character, snch as charge transfer states [159,160]. 

In snbsystem DFT, using the idea of LR allows one to obtain the response 
of a fnlly interacting system (y) by calculating the response of each snbsystem 
{xj) following simple eqnation [152]: 


/ 

where by xl we mean a snbsystem response fnnction that inclndes the in- 
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teraction (coupling) with all other subsystems. The couplings, which are 
manifested by changes in the densities of subsystems induced by the changes 
in the densities of other subsystem, hold valuable physical insight on the na¬ 
ture of the subsystems and the embedding, unaccessible from a supersystem 
calculation. 

In the framework of subsystem TDDFT, the effective potential acting on 
each subsystem is given by: 


= 5vappi{Y,t) + (40) 

and in the frequency domain for subsystems, Eq.(32) transforms to: 

5pi{r,u) = J xj{^,r',u})Svappi{r',u)dr' 

= j X°i{r,r',uj)6viff{r',uj)dr' (41) 

where X/(r, r', u) is the fully coupled response function of subsystem I. The 
term is a measure of the coupling of the given subsystem with 

itself and with the environment induced by the perturbation. Hence, it is 
related to the density fluctuations in all subsystems. Namely, 

Ns 

J 

where, 

Kijir, r, u) = + /xc(r, r, ta) /T(r, r, u) - /r(r, r, uj)5ij (43) 

where /a;c(r, r', ca) is the exchange correlation kernel and /T(r,r',a;) and 


iF/j(r, r, oj ) 5 pj { y \ a;)dr' 


(42) 
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/^(r, r', cj) are the kinetic energy kernels of the full system and for each sub¬ 
system I, respectively. Eq. (42-43) are due to Neugebauer [153]. The kinetic 
kernels, when expressed in the time-regime, have the following structure: 


/r(r,r',t -t') = 




Spirit) 5 pir'^t') 

5‘^T,[pi] 


(44) 


(45) 


5pi{Y,t)5pi{Y',t') 

Using Eq. (42) in the expression for density changes, [Eq. (32)] and trans¬ 
forming to the frequency domain, we arrive at: 


5pi{v,u) = / x?(r,r',a;)(5nappz(r',a;)dr'-Fx?(r,r',a;) ^it'/j(r',r",a;)(5pj(r",a;)dr'dr" 


(46) 

This equation shows how the induced potential couples all density responses 
of the subsystems with each other. The coupling in a given subsystem with 
itself, dehned as xl is readily obtained as [161]: 


(x'l) '(r.r'.w) = (x?) ‘(r.r'.w) -/0,(r,r',w) 


(47) 


Inserting Eq.(47) in Eq.(46) we get: 


(Xi) \r,r',u;)6pj(r',u;)dr' = dvappi(r,u;) + / yKjj(r,r',u;)6pj(r',u;)dr' 


J^i 


(48) 

If the changes in subsystem density are expressed in terms of the coupled 
subsystem response function and the applied potential, governed by Eq.(41), 
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one arrives at a Dyson-type equation given by: 


Ns 

= X'I + Y1 (49) 

where, if the following expression for x'l is used: 

X^ = xUx"iKnX^ (50) 

we obtain: 

Ns 

X/ = X? + X? KijXj (51) 

j 

Eqs. (49)-(51) connect the total, uncoupled and KS responses for the subsys¬ 
tem. It is important to point out that in the derivation of these equations, 
no further approximations where employed and, in principle, the subsystem 
LR-TDDFT method is applicable to any system, provided the ground state 
densities can be correctly described using subsystem DFT with the presently 
available NAKF functionals. 

4.2.1 Applications using subsystem LR-TDDFT 

The original subsystem LR-TDDFT introduced by Casida and Wesolowski 
[152] offered an approximated but efficient way to include the effect of embed¬ 
ding into the the calculation of frequency dependent properties but calculat¬ 
ing the uncoupled response [Fq. (50)] of each subsystem, which includes the 
couplings of the excitations within the subsystem but not with its surround¬ 
ings. When applied to systems where the molecular orbitals are localized 
on the same molecule, this proved to be a very good approximation: for the 
calculation of excitation energies for guanine-cytosine (G-C) and adenine- 
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thymine (A-T) base pairs in the Watson-Crick arrangements, the uncoupled 
version of subsystem LR-TDDFT was able to reproduce excitation energies 
differing by less than 0.05 eV compared to supramolecular calculations. [162]. 
Additionally, the authors were able to extract interaction induced shifts of 
the monomers compared to the isolated molecules, which ranged between 
—0.17 eV and 0.32 eV. It is interesting to note that the presence of hydrogen 
bonds in the dimers did not have an effect on the accuracy of the method: the 
difference between and x'^ is not so much affected by the presence of bonds, 
but by a direct entanglement of the excitations between the subsystems. 

The uncoupled version of subsystem LR-TDDFT, from now on labeled 
FDEu, proved to be very useful for the evaluation of solvatochromic effects, 
where the presence of a solvent results in a shift of an absorption or emission 
band. [163] When explicit interactions between the solvent and solute exist, 
for example in the form of hydrogen bond interactions, such effects are known 
to be problematic to model using implicit models, which only include the 
dielectric medium effects [164]. Using FDEu one can apply even a further 
approximation, by restricting the orbital space to the relatively smaller solute 
system, while including the solvent only as a frozen density for both the 
ground state and time dependent FDE calculation. In such cases FDEu also 
offers the advantage of only supplying the excitations of the subsystem in 
interest, where in supramolecular calculations the calculation cost is driven 
further up by the large number of excited states that needs to be included. 
For further reduction of the computational cost, one can generate the density 
of the solvent as a superposition of spherically symmetric atomic charges, 
though it has been shown that this method can lead to large errors when the 
molecules within the environment interact strongly by, for example, forming 
hydrogen-bonded chains. [40] In such cases it is recommended to hrst perform 
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a supermolecular calculation at a lower level of theory. An application on 
the calculation of excitations of an acetone molecule solvated in water, it was 
found that FDEu was able to produce the correct trend for the energies of the 
valence excitations with increasing cluster size and the solvatochromic shift 
was found to be in excellent agreement with experiment. [165] An extensive 
analysis of the structural effect due to the dynamics in solution and the 
electronic effect due to the frozen densities was performed on a study of 
aminocoumarin C151 solvated in water and n-hexane [87] showed that the 
electronic effect is more important and a relaxation of the subsystem densities 
using freeze-and-thaw cycle can be important. In particular, polarization 
of the density becomes important in cases where there is direct hydrogen 
bonding between the frozen and embedded subsystems, as in case of water. 

While FDEu proved to be successful on many occasions, it also raised 
the question of when and start to diverge sufficiently for FDEu to 
fail. [166] Indeed, in some cases the couplings between the subsystem re¬ 
sponses can become dominant. An evident example is a pair of coupled 
chromophores, where excitation energy transfer occurs. If the two chro- 
mophores are identical, one can extract information about the excitation 
rate transfer from the transition-dipole-transition-dipole interaction, which 
expresses itself in the optical absorption spectrum as a splitting of the cou¬ 
pled excitation energy. For this purpose, the subsystem TDDFT version 
was extended [153] to include the calculation of between selected inter¬ 
subsystem excitations. Even when all couplings are taken into account, the 
dimension of the total problem is still smaller than in a supermolecular cal¬ 
culation since there are no orbitals dehned for the supersystem and hence no 
inter-subsystem orbital transitions. However, the computational advantage 
of FDE is taken further into acconnt by introducing further approximations. 
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One first solves the FDEu problem, which will contain most of the informa¬ 
tion on the environmental effects due to the embedding potential. One can 
proceed by only calculating ~x^ for the exciton-like couplings between local 
transitions. [153,167,168] For the case of benzaldehyde dimer [153], for inter- 
molecular distances larger than 5 A, where the excitonic coupling is small, 
both the FDEu and FDEc methods produced similar results. For the shorter 
distance of 4 A, however, a difference of 0.15 eV was observed. It should 
be noted that FDEc does not include couplings to excitations with partial 
charge-transfer characters, as the transition densities are expressed in terms 
of the monomer orbitals. The method can be extended to solvent effects, 
including selected solvent excited states. [38,73] The selection of which exci¬ 
tation should be included in the coupling of the FDEc method, the authors 
used the sum-over-states (SOS) expression of the polarizability. The selection 
was made by ranking the excited states of the solvents in decreasing order of 
contribution to the SOS polarizability and choosing the first k states for re¬ 
producing the full polarizability until a threshold of for the desired accuracy 
was achieved. Additional properties for which subsystem LR-TDDFT can 
be applied include polarizabilities and optical rotation parameters which can 
be used for the calculation of oscillator and rotatory strengths [169], induced 
circular dichroism, [163,166] vibronic spectroscopy. [163,170,171] and reso¬ 
nance Raman spectra. [172] We refer the readers to Ref. [25] for a detailed 
overview of specihcs of the methods and applications. 

A recent development is employing the linear response formalism to per¬ 
form a WF-in-DFT embedding, as opposed to DFT-in-DFT embedding as 
discussed until now, in order to model local electronic excitations. [13,104, 
173-176] The formal foundations to make this possible were laid in Ref. [177], 
where the expression for the local embedding potential as a functional of 


50 



charge densities was derived in a general case of embedded interacting Hamil¬ 
tonian. An excellent review on the subject can be found in Ref. [13]. 

4.3 Real-time subsystem TD-DFT approach 

Real-time TDDFT (rt-TDDFT) involves the direct integration of the TDKS 
equations [Eq. (30)] in time. The main advantages of this approach are that it 
allows to model the response of the system beyond the linear response, obtain 
optical absorption spectra without the use of virtual Kohn-Sham orbitals 
and study transport properties in real time. It can also be extended to full 
dynamics of the system by also evolving the nuclei in time and for very large 
systems it even has computational advantages as it scales as 0{N^). A full 
discussion of the rt-TDDFT method lies beyond the scope of this review and 
we will concentrate on the most popular approaches for solving Eq. (30) and 
how they can be extended to various flavors of subsystem rt-TDDFT. The 
rt-TDKS equations are solved in the following steps: 

1. A ground state KS-DFT calculation is performed, obtaining the KS 
orbitals 4>j{r). These orbitals constitute the initial condition for the 
time evolution of the system. In principle, the initial condition is not 
required to be a ground state, as long as the initial state of the nonin¬ 
teracting system <Fo can be represented as a single Slater determinant. 
For instance, one can construct an initial an excited state by promoting 
one of the electron to a virtual orbital and performing a constrained 
DFT calculation with the new set of occupied orbitals. 

2. The time dependent effective potential is set up as 

ns[p](r,t) = nappi(r,t)next(r,t) j + v^c[p]{r,t) (52) 


51 



where we adopt the adiabatic approximation for the XC functional, 
which allows us to use the static XC functional from the ground state 
calculation. The choice of the time dependent applied held depends on 
the nature of the application and will be discussed below. 

3. The integration in time of Eq. (30) is performed, usually using the 
Crank-Nicolson algorithm [178], by solving the linear set of equations 

([i + ~ ^ 

(53) 

The At is needs to be taken very small, around 1-2 attoseconds, to en¬ 
sure numerical stability. The advantage of the Crank-Nicolson scheme 
is the preservation of the energy and the orthonormality of the KS 
orbitals. 

Currently there are two havors of subsystem real time TDDFT implemen¬ 
tations, with exact embedding theories [14,15] and with FDE [154]. The hrst 
implementation of fragment-based TD-PDFT was introduced by Mosequera 
et al. [14], where a formal proof of the unique mapping between the total 
time dependent density and the time dependent partition potential is given. 
As in the ground state case, the total system is divided into fragments with 
the restraint that the sum of the fragment densities at each point of space 
and time equals the total density 

p{r,t) = ^pi{r,t). (54) 

I 

One starts by first performing a ground state PDFT calculation in order 
to obtain the initial fragment KS orbitals. The fragment KS orbitals are 
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propagated in time using the following time-dependent fragment equations 


f)]0^(r, t) (55) 

where the effective fragment potential Vs^i{r,t) contains the partition poten¬ 
tial np(r, t). The partition potential is determined at each time step from the 
total density of the system: 

V-|p(r,«)V'tip(r,()| = ^ ** +^(V-Q.,jI'iipKr, i)-V- \p,[i]^] (r, t)V-ii.j[-tip](r, i)]) 

(56) 

where Q is dehned using the density matrix T and the current density j 

Qs,i =-itr{Tsj{t)[j{r),f]} (57) 


The calculation of the time dependent partition potential at each point of 
space thus requires an independent calculation of the total density. For this 
reason, the total density is propagated simultaneously with the densities of 
the fragments. 

In the TD-PFET method developed by Huang et ah [15], the necessity 
to calculate the total density at each time step is omitted, as they derived a 
formal equation for the time dependent embedding potential from the action 
formalism 


S^tot ["i^emb] 
(5nemb(r,t) 



['^emb] 


(5'htoi,r['yemb] \ 
(5nemb(r,t) / 


(58) 


When the embedding calculations are carried out in such a way that all 
subsystems are treated within TD-DFT using adiabatic local density ap- 


53 



proximation (ALDA), this formal expression can be approximated as: 

^eff[Ptot] A '^xc[Ptoi] ^niPioi] Aon,tot Ad (^9) 

Optot 

Given that at Wemb equals zero, this allows to calculate the embedding 
potential iteratively 

t) = t) - (60) 

OPtot 

where n is the iteration number. To perform a TD-PFET calculation, one 
starts with a ground state PFET calculation to obtain the fragment densities 
and Kohn-Sham orbitals, as well as Vemb(i‘), which is used as the initial guess 
for Vemb(i‘G) at to- The fragment Kohn-Sham orbitals are propagated in 
time using the Crank-Nicolson method with a trial Vemb(i’G) and used to 
construct the fragment and total densities at t -|- At. The effective potential 
and total Kohn-Sham orbitals associated with the total electron density are 
then obtained using the time-dependent version of the Zhao-Morrison-Parr 
method [179,180], which involves the solving of a set of coupled equations 
self-consistently. From these, one can obtain and a new trial nemb(i'G)- 
This process is repeated iteratively at each time step until conversion of the 
total electron density. In other words, omitting the necessity to perform a 
supramolecular rt-TDDFT calculation as done in TD-PDFT comes at the 
cost of two nested iterative calculations for obtaining the time dependent 
density and the effective potential. 

An alternative method was introduced by Krishtal et al. [154], where the 
subsystem DFT formulation is extended into a real-time subsystem TDDFT 
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by evolving the TDKS equations in time for each subsystem 


--V^ +viir, tM 




( 61 ) 


The time dependent effective potential of each subsystem is dehned as 




5J[p] , 5E,,[p] , 5t[p] 


+ 


+ 


6p{r,t) 6p{r,t) 6p{r,t) 


St[pi] 

6pi{r,t) 


(62) 


and is equal at t = to to the (r) in the solution of the ground state subsys¬ 
tem DFT. Only Vext(i', t) contains an explicit time dependence while the other 
terms only depend on time through the density. If all subsystems are propa¬ 
gated simultaneously and the total density is updated at every time step, the 
full rt-TDKS solution is recovered within the accuracy of the NAKE. If the 
applied external held is sufficiently small for the density response not to di¬ 
verge from the linear response regime, one can draw a parallel to the response 
functions from linear response TDDFT, discussed in more detail in Sec. 4.2. 
In that case, the solution of each subsystem includes the full coupled response 
of the subsystem xj [161], as will be discussed below in Section 4.2. If Eq. 
(61) is only propagated for one of the subsystems, while keeping the other 
subsystems frozen, the solution will only include the uncoupled subsystem 
response xl- The subsystem rt-TDDFT method was implemented into the 
Quantum Espresso package using plane wave basis sets, as described in 
Section 2.2.2. The solution of the equations is performed in a similar fash¬ 
ion to the rt-TDDFT, by applying the 3 steps described at the beginning of 
this Section to each subsystem. For the case of the full coupled calculation, 
the total time dependent density is reconstructed from the densities of the 
subsystems at the end of step 3. 
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The three subsystem rt-TDDFT methods are summarized in Table 2 
w.r.t. main quantities, the required input at the initial time to, the required 
input at each time step and the output quantities. 

4.3.1 Electronic spectra 

Electronic spectra can be obtained with rt-TDDFT by applying a short laser 
pulse of strength e to the Kohn-Sham orbitals of the system, propagating 
them in time and Fourier transforming the time-dependent dipole moment 
fi{t) into the frequency domain. The pulse can be either applied by shifting 
the Kohn-Sham states [181] 


0i(r,t = 0+) = e*^''0i(r,t = 0 ) (63) 

where e is the field strength or by adding an explicit time dependent potential 
in the form of a very narrow gaussian in time that integrates to e. The 
second option is not feasible in periodic calculations but both alternatives 
produce identical results for molecular systems. If e is sufficiently small, the 
density response will be conhned to the linear response regime and the results 
will be directly comparable to those obtained the using the linear response 
TDDFT formalism, discussed in Section 4.2. By applying a stronger held, 
one can straightforwardly study effects beyond the linear response, which 
is one of the main advantages of rt-TDDFT. This, however, requires extra 
care in the numerical time integration of Eqs. (61), such as the use of very 
small time steps (< 1 as) and a predictor-corrector scheme for the Crank- 
Nicholson algorithm. All applications reported here were performed using a 
small electric held. 
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The oscillator strength at each frequency uj is given by 

Sk{^) = - — / sin(a;t)e"^*[/ifc(t) -(64) 

efcvr J 

where 7 is a small damping factor associated with the rj factor in Eq. (35). 
The resolution Aca of the spectrum depends on the simulation time and 
usually a simulation of at least 10 fs is necessary. At the present, no results 
on electronic spectra have been reported using the TD-PFET and TD-PDFT 
methods, so we will concentrate on results obtained using subsystem rt- 
TDDFT method [154]. 

An example of what kind of information that can be obtained from a 
subsystem rt-TDDFT method is a dimer of stacked benzene and fulvene 
molecules, separated by a distance of 5 A. The optical spectrum is obtained 
by applying a laser pulse of 0.01 Ryd A and evolving the Kohn-Sham orbitals 
of each of the subsystems in time for 8000 steps with a time step of 2 as. The 
dimer was placed in a supercell of 31.0 x 32.5 x 37.8 a.u.^ to avoid periodic 
interactions. The calculation was performed with the PBE functional [182], 
the LC94 [183] functional for NAKE and ultrasoft pseudopotentials from the 
GBRV [184] library with a kinetic energy cutoff of 55.0 Ry and density cutoff 
of 660.0 Ry. Fig. 9 compares the results obtained using rt-TDDFT and sub¬ 
system rt-TDDFT. The spectrum of the dimer in the subsystem rt-TDDFT 
method is obtained by simply adding the spectra of the two subsystems, 
since the time dependent dipole vector is additive across the subsystems. 
Since benzene and fulvene do not absorb in the visual region, the spectrum 
is shown for the 3.26-9 eV range. Higher energy frequencies are not con¬ 
sidered here. FDE succeeds in reproducing the TDDFT results with slight 
differences resulting from the approximate kinetic energy functionals used in 
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the embedding potential. All excitations under 8 eV are reproduced within 
the 0.01 eV accuracy, with only slight deviations in the oscillation strength. 
At higher energies region, the deviations are slightly larger: this could also 
be an attribute of the numerical accuracy of the rt-TDDFT implementation, 
since higher frequencies are reproduced at shorter simulation times and must 
be commensurate to the time-step. 

[Figure 9 about here.] 

While reproducing the supersystem result is an important feature of any 
subsystem method, the real interest lies in comparing the optical spectra of 
the subsystems with the optical spectra of the isolated entities. As can be 
seen from Fig. 10, where the optical spectra of the isolated fulvene molecule 
and the fulvene molecule in the dimer are depicted, the embedding potential 
generated by the presence of the benzene molecule results in a shift to lower 
energies. The shift is present for all excitations, but is especially pronounced 
for the excitation lightly higher than 6 eV. Furthermore, the embedding 
potential generated by the presence of the benzene molecule has a clear effect 
on the intensity of the peaks, which is especially pronounced for the excitation 
at 6.44 eV in the isolated fulvene spectrum. 

[Figure 10 about here.] 

An important question is the extent to which the full dynamic response is 
vital for the reproduction of the interaction induced shifts and changes in 
the oscillation strength. This question has been broadly discussed for the 
subsystem formulation of linear response TDDFT. In the hrst formulation 
of the method, one performs a linear response TDDFT calculation on one 
subsystem while statically including the density of the other subsystem(s) in 
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the embedding potential but restricting the response to the active subsystem 
only. [152,162], Several studies have shown that this approach is sufficient 
in reproducing optical spectra [87,165], Raman spectra [172] induced circu¬ 
lar dichroism [163,166] and electron-spin-resonance hyperhne couplings [185] 
even in the presence of hydrogen bonding, as long as there are no explicit cou¬ 
plings in the excitations between the systems. Later on, Neugebauer [153] 
introduced an approach to include couplings between selected excitations. 
This effect is clearly seen when comparing the optical spectra of the embed¬ 
ded fulvene molecule from a coupled and uncoupled calculations, as depicted 
in Fig. 11. In the coupled calculation, both molecules were subjected to a 
pulse and integrated in time, with a total density update after each time 
step. In the uncoupled calculation, the density of the benzene molecule is 
kept frozen during the full length of the simulation, which is equivalent to 
performing the uncoupled version of the linear response FDE calculation. 

[Figure 11 about here.] 

One can see that both methods produce very similar results for most excita¬ 
tion energies with the exception of the excitation energy at 6.44 eV, where 
there is a difference in both the interaction induced shift as the associated 
oscillator strength. The reason for this difference is that the benzene and 
fulvene molecules are strongly coupled at this frequency [154], While the un¬ 
coupled calculation succeeds in reproducing the results for most excitations 
from an electronic effect through the presence of the density of the other sub¬ 
system in the embedding potential, it underestimates the effect generated by 
the other subsystems for the coupled excitations. In this particular case, 
the interaction induced shift is not as strongly pronounced in the uncoupled 
calculation as in the coupled calculation and resembles more the excitation 
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energy in the isolated fulvene molecule. A similar effect is found in the dimer 
with a shorter intermolecular distance of 4 A, where an interaction induced 
shift was observed in the coupled calculation for the excitation energy at 
4.99 eV in the fulvene molecule and 6.78 eV in the benzene molecule. This 
interaction induced shift in the uncoupled calculation amounted only to half 
of the value. It is also interesting to note that the coupled calculations pro¬ 
duce lower oscillation strengths for the coupled excitations, compared to the 
uncoupled calculations and the isolated molecules. This is compensated by 
higher oscillation strength values at the uncoupled excitations, in order to 
satisfy the sum rule of the spectrum. 

4.3.2 Transport properties 

An important advantage of rt-TDDFT is the ability to study processes in 
real time, where the system evolves due to an outside perturbation such as 
an applied field or an excitation, i.e. thinks of processes such as electronic 
conductance, excitation transfer and charge transfer. Since subsystem rt- 
TDDFT has the restriction of constant charge on the subsystems, charge 
transfer cannot be studied straightforwardly, although it is possible by using 
the restriction to our advantage and calculating electronic couplings between 
charge-localized, diabatic states [128]. The real limitation for the study of 
charge transfer is the lack of good quality kinetic energy that can adequately 
describe overlapping densities and XC functionals that can reliably describe 
charge separated states. With the exact kinetic energy functional, subsys¬ 
tem DFT would reproduce the total DFT result and therefore also a charge 
transfer reaction. However, this would inevitably result in a delocalization of 
the subsystem electron densities and while the total electron density would 
reproduce the KS-DFT one, unless self-interaction corrected XC functionals 
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are used it would likely be incorrect. 

To demonstrate the ability of subsystem rt-TDDFT, we will illustrate 
the explicit transfer of excitation energy between two coupled systems as 
discussed by Krishtal et ah [154], An Na 4 cluster, consisting of two Na 2 
molecules, is a well studied system [15,186,187] which couples strongly at 
the excitation energy of 2.18 eV. For this purpose, two Na 2 molecules were 
placed at a distance of 6.6 A. At the start of the simulation, an electric field 
in a direction along the Na-Na bond is applied to only one of the subsystems, 
noted as the donor, with a frequency corresponding to an excitation energy of 
the sodium dimer of 2.18 eV. Both of the subsystems are evolved simultaneous 
for 100 fs with a time step of 2 as. The full cluster is placed in a supercell 
of 22.7 X 43.5 x 22.7 a.u.^ to avoid periodic interactions. The calculation 
was performed with the PBE functional [182], the LC94 [183] functional for 
NAKE and ultrasoft pseudopotentials from the GBRV [184] library with a 
kinetic energy cutoff of 55.0 Ry and density cutoff of 660.0 Ry. 

Fig. 12 depicts the evolution of the dipole moment along in the direction 
of the applied held. As one can see, at the beginning of the simulation, the 
dipole moment of the donor reacts to the applied held and the donor molecule 
is excited. The system proceeds in transferring the excitation energy in a 
periodic fashion between the donor and the acceptor molecules, as can be 
clearly seen in Fig. 12. Note that the acceptor subsystem has no direct 
applied held in its potential and only experiences the excitation due the 
embedding potential. The alternate beating of the dipole moment of the 
two subsystems as a direct consequence of excitation energy transfer and 
a full cycle, between the maxima of two consequent beatings on the same 
subsystem, is the rate of the excitation transfer. For the particular separation 
distance of 6.6 A, it equals to 19.1 fs. The excitation energy transfer rate is a 
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well studied phenomenon, related to the Forster and Dexter energy transfer 
theories. Forster [188,189] described excitation energy transfer in a simplihed 
model using a dipole-dipole interaction of the transition densities of the two 
systems, with an ^ dependence. This model is valid in the long range 
where the Coulomb interaction is dominant. For shorter range, Dexter [190] 
extended the model to include higher multipole order and exchange effects. In 
linear response formalism, excitation energy transfer is studied by calculating 
the excitation energy transfer couplings explicitly. Needless to say, in order 
to model this effect using subsystem TDDFT both in the linear response and 
real-time formalisms, a coupled formulation of the method is needed. [73,153]. 

[Figure 12 about here.] 


5 Subsystem prospective of many-body in¬ 
teractions 

The calculation of such long-range interactions as van der Waals in molecules 
pose a major challenge. This is due to the failure of the Local Density Ap¬ 
proximation (LDA), which makes use of a homogeneous electron gas model, 
and the Generalized Gradient Approximation (GGA), which is semilocal in 
nature. Because of this, KS-DFT calculations carried out with semilocal XG 
functionals is only partially able to account for interactions that are nonlo¬ 
cal in nature, such as dispersion interactions [191,192] and all long-ranged 
interactions originating from the correlated part of the energy functional. 
However, in those cases where the electron density is nonvanishing in the 
region separating molecular fragments, semilocal approximations can still ef¬ 
fectively account for long-ranged interactions due to the correlation energy 
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(see for example Ref. [193]). 

In the following, an in-depth analysis of existing [194] and novel [115] 
subsystem formulations of the adiabatic-connection fluctuation-dissipation 
theorem (AC-FDT) in DFT will be described with a particular focus on how 
these formulations handle density overlap between the subsystems. Practical 
implementations of these methods as well as related ones [195-197] will be 
discussed and preliminary results of the van der Waals inclusive subsystem 
DFT method will be presented [115]. 

In the AC-FDT, the correlation energy is related to the response functions 
of the fully-interacting system, y, and the one of the noninteracting KS 
system, xo- The key is scaling the electron-electron interaction by a coupling 
constant. A, in a way that when A = 0 one recovers the KS system and when 
A = 1 the interacting system is recovered [198,199]. Hence, if one starts from 
the ground state energy equation for a particular value of A, given by: 

= (I'SlJ/t'to) (65) 

it is clear that Eq = Eks for A = 0, i.e., the Kohn-Sham energy and Eq = Eq 
for A = 1, i.e., the exact ground state energy of the system. Thus we get: 

Eq = Eks + j 

Through the fluctuation-dissipation theorem [200], Eq. (66) is related to the 
response functions yielding the following expression for the correlation energy 

Ecorr = - dAdrdr'^- - —Im[x^{r,r,u)-xo{r,r,u)]. (67) 

According to the prescribed scaling of the electron-electron interaction, 
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is given by the Dyson equation: 


X''(r,r',n)) = Xo(r,r',w)+ 





( 68 ) 

Thus, van der Waals interactions can be estimated using the AC-FDT in an 
effective manner [201]. 


5.1 Adiabatic-connection fluctuation-dissipation theo¬ 
rem for subsystems with nonoverlapping densities 

The Hamiltonian for a model system comprising of two subsystems, A and 
B, could be separated into three components: 


II = IlA + ffB + VaB 


(69) 


where IIa and Hb are the components of isolated subsystems and Vab is the 
interaction which depends on the distance between them and is given by: 


ri-ral ri^ Y2a 

Herein, ri /2 belongs to electrons in subsystem 1/2, and Yix is the electron- 
nuclear separation between electrons of subsystem i and nuclei X. The 
Coulomb interaction can be scaled by an adiabatic-connection-like param¬ 
eter, in the same spirit as AC-FDT. Upon carrying out steps similar to 
Eq. (65-68), Dobson & Gould [194] derived the nonretarded Lifshitz formula 
which is closely related to the generalized Casimir-Polder equations for two 
distinguishable systems (which was derived already with second-order per- 



Vab{^i 
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turbation theory [197,202]), 


^+oo 

-Edisp = - 7 ^ dca 

27r Jq 


Xi{r'i,ri,iu)x2{r2,r'2,iuj) 
|ri - r2||r'i - r'2| 


dridr2dr'idr'2. 


(71) 


Dobson & Gould made two assumptions: (1) the use of response functions 
from isolated monomers in the full-potential approximation (i.e. the response 
functions do not depend on the AC-FDT coupling parameter), and (2) that 
the subsystem’s response function is given by (in a simplihed notation) 


X = Xi + X2 + / dridr2. (72) 

J |ri-r2| 

The above equation is an approximation because it accounts for the inter¬ 
actions between the response functions of the two subsystems in the RPA 
approximation, rather than with the full interaction of Eq. (43). 


5.2 AC-FDT for density-overlapping subsystems 

In the subsystem DFT scenario, the effective interaction between the sub¬ 
system response functions differs from the RPA approximation in the fact 
that it includes XC and KE terms. This is a pivotal point differentiating the 
approximate Eq. (71) form the exact one. For example, the so-called overlap 
effects are not included in the dispersion energy in Eq. (71) [202]. Neglecting 
these effects is at the origin of the inclusion of damping factors in common 
van der Waals corrections to semilocal KS-DFT [203]. 

An exact expression for the correlation energy of the interaction between 
a collection of subsystems can be derived form the subsystem DFT dehnition 
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of additive and nonadditive correlation energy functional. Namely, 


Ns 


Ec — [pi] + [p/, P/7, , PtVs] 


(73) 


Long Range 


Short Range 


As the correlation energy is related to the response functions by the AC-FDT, 
and because of Eq. (39), we have that 


Axa = x" - Xo = Axf+ Ayr. (74) 


When the system is separated into smaller subsystems, the coupled response 
function for each subsystem from Eq. (49)] and the uncoupled response 
function for each subsystem [y/’^, from Eq. (50)] must also depend on the 
value of the coupling constant A, and the difference between these two re¬ 
sponse functions will, in turn, be related to the correlation effects in each 
subsystem. The correlation contributions from each subsystem should then 
sum up to give the total correlation of the entire system, which, as has been 
mentioned already, resides in the nonadditive component of the A dependent 
response functions giving rise to the following equation: 

Ns 

Axr = E (xh - xp) 

I 

Ns 

= ( 75 ) 

/ 

Hence, using the expression for correlation energy, given in Eq.(67) above. 
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and Eq.(49), we arrive at the following equation: 


^nad _-J 

" 2n 



dxidx; 


/o ir 12% xf\‘^)KhHxf{‘^) 


ri - T2\ 


(76) 

The above equation achieves the goal of expressing the exact correlation 
energy of interaction from subsystem quantities (i.e. the subsystem response 
functions) and appropriate kernels of interaction: Coulomb, XC and KE 
kernels. It generalizes the result of Dobson & Gould for nonoverlapping 
subsystems in a soft way, i.e. Eq. (76) still resembles the generalized Casimir- 
Polder formula, Eq. (71). 

As all exact equations, it is of impossible evaluation. Thus approxima¬ 
tions need to be introduced, and will be discussed in the following section. 


5.3 Practical implementation 

The nonadditive correlation energy determined with Eq. (76) is still compu¬ 
tationally prohibitive for realistic systems, as the coupled response functions 
must be obtained solving the Dyson-type equation in Eq. (49) which in the¬ 
ory requires operations, where N here is the size of the entire supersys¬ 
tem [161,167]. One approximation is achieved by considering a perturbative 
solution to the Dyson-type equation Eq. (49) as 


Ns Ns 



Jy^I Jy^I, K^J 


The above expansion shows how the coupled subsystem response functions 
contain terms of polarization by other subsystems, which come at the n- 
body level. Retaining only the first term of the expansion (e.g. a pair-wise 



approximation) Eq. (76) is simplified. The nonadditive correlation energy is 
now expressed in terms of x'l only - i.e. the response polarization is neglected, 
and only the polarization arising in the ground state is accounted for through 
the dependence of x'l fo the subsystem orbitals. The next approximation 
stems from the fact that Kjj is not known exactly. Motivated by the success 
of RPA for noncovalent interactions, we approximate Kjj by neglecting the 
frequency-dependent exchange-correlation kernel as well as the kinetic energy 
contributions to the kernel (/xc, /t, and /rf) 


K 


ij 


K 


A,RPA 
IJ 



(78) 


Since the RPA kernel is frequency independent, it is computationally much 
more efficient to solve. Finally, the last approximation is the full-potential 
approximation for the response functions. We summarize the various ap¬ 
proximations and their acronyms as proposed in Ref. [115] in Table 3. 


[Table 1 about here. 


For only two interacting subsystems {Ns = 2), application of the GCP)* 
approximation and subsequent integration over A (which gives a factor 1) 
yields 


^nadd _- 

27r 



/ dxdx' dcux)‘(xi, X, u)t^^X2{^', ^ 2 , oj) 


dxidx^ 


|r-r'|' 


Ti - r2 


(79) 

Clearly, the approximations that led us to Eq. (79) are similar to the ones 
leading to Eq. (71). Thus, the two equations resemble each other with one key 
difference, x'l ¥" Xi- With that, we should expect Eq. (79) to yield somewhat 
superior dispersion energies than Eq. (71), provided that the ground state 
density and energy are satisfactorily characterized by the FDE calculation. 
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When implementing Eq. (79), for closed-shell systems, the spectral repre¬ 
sentation of the response functions, assuming real solutions of the TD-DFT 
eigenvectors, is given by [204] 


Xy(r,r',a;) = ^ 


Aul 




E + y:ux: + y:),^ 

(ia)i,{jb)i 

0*(r)0a(r)0j(r')0b(r'), (80) 


where 0j(r) are the KS orbitals of subsystem I. Subscripts denote 

occupied orbitals and a, b, c, d virtual orbitals. (X“ -1- Y^)ia is the projection 
of the sum of the excitation (X) and de-excitation (E) TD-DFT eigenvector 
for the n-th excited state of subsystem I. The associated eigenvalue is given 
by The Hartree-XC kernel Ku used to determine xj is given by Eq. (43) 
with I = J [152,165]. Using Eq. (80) and employing the GCP“ approximation 
(see Table 3) we can express the difference of the coupled and uncoupled 
response functions accordingly, and the AC-FDT formula takes the form 

Ns Ns ojvi ojvj 
I J>1 {n)i {m)j 

( A '„"+ y:ux: + y:)„{xi ++ 1 -),, ( si ) 

{ia\ld) {jh\kc). 

The above formula scales as x iVj for each pair of subsystems considered. 

A pilot implementation of Eq. (81) has been carried out in a local version 
of the ADF computer package [205], and was presented in Ref. [115]. The ap¬ 
plication to a subset of the S22 set gave encouraging results, e.g. the binding 
energies of molecular dyads were improved signihcantly over FDE computed 
with semilocal nonadditive correlation functionals. 


« + 


OJVI OjVl OjVj OjVj 

EEEE 

{ta)i {jb)i {kc)j {ld)j 
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6 Future directions 


The future of density embedding methods holds promise. This is because 
partitioning the supersystem into subsystems can be achieved by many dif¬ 
ferent approaches. Thus, the method can capitalize on a plethora of avenues 
of research. 

In Section 2.2.2, we have advocated the possibility of using FDE for the 
description of periodic systems. Recently [144], the fact that FDE does not 
impose orthogonality between the subsystem orbitals and yields localized 
densities was exploited to calculate dipole moments of periodic systems us¬ 
ing the machinery proper of molecular systems. This provides us with an 
additional motivation to pursue an ongoing project in our lab which regards 
devising a subsystem-specihc hrst Brillouin zone sampling (the so-called k- 
point sampling). This shows that FDE is a very flexible method (e.g. the sub¬ 
systems can be treated at different level of theory). However, when ground 
states are considered, there are no doubts that in order to improve the ap¬ 
plicability of FDE, new and more accurate nonadditive Kinetic Energy func¬ 
tionals need to be formulated - a strong overlap between subsystems leads 
to failure of FDE when semilocal NAKE functionals are employed. Efforts 
in this direction began with the introduction of nonlocal KE functionals [33] 
and are still ongoing by several groups [39,99,100,206-209]. Despite this, a 
true breakthrough has still to come for the noninteracting KE functionals. 

When departing from the concept of density embedding, the Pauli Block¬ 
ade (PB) method (reviewed in Section 1.3) provides an elegant solution to the 
strong-overlap failure of approximate NAKEs. This method, especially when 
coupled with a correlated wavefunction treatment of the subsystems [65-67], 
has already found important applications [210]. Other types of embedding 
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include the so-called density-matrix embedding methods [211-214] having 
the advantage of providing an exact embedding framework for uncorrelated 
(Hartree-Fock, and Kohn-Sham) as well as correlated wavefunctions. The 
future holds the possibility of using PB for calculating excited states (and 
their properties), as the linear-response theory associated to the PB method 
has not been developed yet. 

The subsystem DFT theory of many-body interactions [115] presented in 
Section 5 has been implemented in the RPA approximation and using only 
pair-wise terms in the response function perturbative expansion of Eq. (77). 
In addition, the full-potential approximation was also used (i.e. the response 
functions are assumed independent on the adiabatic-connection parameter). 
These choices are very restrictive, as the two-body approximation is known 
to break down in many instances [215,216], and the RPA approximation is 
responsible for overestimated correlation energies [115,217,218]. To go be¬ 
yond the RPA approximation, nonlocal XC and KE functionals are needed 
[115,219-222] carrying with them an apparent issue of computational com¬ 
plexity due to the inherent double-integration. Thus, more work is needed 
to fold in nonlocality in a computationally feasible way, such as by using ap¬ 
propriate htting techniques [223], ad-hoc corrections [224], or approximating 
the nonlocal kernel with computationally feasible ones [220,221]. Lifting the 
full-potential approximation to Eq. (76) will also be explored and its effect 
on the calculated nonadditive correlation energies will be assessed. 

One natural question is whether the method arising from Eq. (79) coupled 
with an FDE calculation of the ground state is capable of reproducing po¬ 
tential energy surfaces (PESs) and not just single point calculations [115]. It 
is known that certain functionals [225] can fortuitously reproduce the bind¬ 
ing energies in the S22 set without actually being able to reproduce the full 
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PESs. Because the currently available NAKE functionals are unable to de¬ 
scribe strongly overlapping subsystem densities, we expect this method as 
implemented with GGA NAKE functionals to fail for short intersubsystem 
separations and for those dimers which are already not well characterized in 
their ground states. Work in this direction is undergoing in our group, and 
will also consider the possibility of adding a van der Waals correction to KS- 
DFT compnted with Eq. (79), using the KS-DFT energy as the nncorrected 
energy as opposed to the FDE energy. 

The real-time TD-DFT implementations (from our group for FDE [154], 
as well as from others for other flavors of subsystem DFT [14, 15]) pro¬ 
vide a unique avenue for exploring excited states of embedded systems and 
their properties with the possibility of going beyond the linear-response 
regime [226,227]. Snpplementing the dynamics of the electrons with the 
motion of the nnclei leads directly to the formulation of a snbsystem nona- 
diabatic Ehrenfest dynamics. The applications of the rt-TD-DFT code go 
even beyond this. For example correlation energies of interaction [e.g. form 
Eq. (76)] can be calculated at imaginary freqnency as pioneered by Marques 
et al. [228,229], rather than by precalculating the linear response fnnction 
with the aid of a set of virtual orbitals as customarily done in molecular 
codes [115,196,217,230]. Due to the time-propagation of the snbsystem or¬ 
bitals, a rt-TD-DFT computation of the nonadditive correlation energy in 
Eq. (76), can be done departing from the perturbative pair-wise expansion 
of the response fnnctions. The longer the time propagation, the more the 
snbsystems will conple. This provides us with a very intuitive way to coarse 
grain in the time-domain which is similar in spirit to coarse graining in the 
imaginary freqnency domain [196,230,231]. 

In conclusion, the fnture looks both bright and busy for the research 
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groups involved in density embedding. 
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evaluate the density and the potential.97 

3 (a): SCF procedure used in [82] as alternative to the freeze- 

and-thaw scheme, (b): MPI architecture of a flexible and 
parallel Subsystem DFT implementation.98 


4 The 3-FDE scheme attempts to allow FDE to simulate frag¬ 
ments connected with covalent bonds, such as polipeptides. 

The total electron density is p = p/ J- pn — Pcap- Reprinted 
with permission from ” Density-functional theory approach for 
the quantum chemical treatment of proteins”, C. R. Jacob and 
L. Visscher, J. Chem. Phys. 128, 155102 (2008). Copyright 
2015, AIP Publishing LLC.99 

5 Isosurface plot of the density difference between a supermolec- 
ular KS calculation and a periodic subsystem DFT calculation. 

(a): Methane on Pt(lOO), top conhguration; (b): Water on 
Pt(lll), horizontal conhguration; (c): Water on Pt(lll), ver¬ 
tical conhguration. Reprinted with permission from ” Periodic 
subsystem density-functional theory”, A. Genova, D. Ceresoli 
and M. Pavanello, J. Chem. Phys. 141, 174101 (2014). Copy¬ 


right 2015, AIP Publishing LCC.100 

6 OH* (H20)ii carried out with the novel FDE implementa¬ 
tion. Spin-density isosurface plot (cutoh of 10“^) obtained 
from EDE and regular KS-DFT are compared. FDE produces 
a (correctly) localized spin-density on the OH* radical.101 
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Schematic representation of the hybrid KS-Frozen environ¬ 
ment method used by Kodak et al. to simulate molecular 
dynamics of solvated systems. Reprinted with permission from 
’’Hybrid ab initio Kohn-Sham density functional theory/frozen- 
density orbital-free density functional theory simulation method 
suitable for biological systems”, M. Kodak, WLu, and J. Bern- 
holc, J. Chem. Phys. 128, 014101 (2008)]. Copyright 2015. 

AIP Publishing LLC. 

Oxygen-oxygen g{r) distribution for several combinations of 
kinetic/exchange-correlation functionals as calculated by lan- 
nuzzi et al.. Reprinted from ’’Density functional embedding 
for molecular systems”, M. lannuzzi, B. Kirchner, and J. Hut- 
ter, Chem. Phys. Lett. 421, 16 (2006). Copyright (2015), with 
permission from Elsevier. 
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non-frozen fragments 


frozen fragments 

electron density Is kept frozen 

aveilable options; 

• Include basis functions ? 

' use Rtted or exact density 7 








relaxed frozen fragments 


electron density is updated In 
freeze-and-tbaw cycles 

available options: 

• number of freeze*and-thaw 
cycles 

• xc functional used In 
freezo*and*thaw cycles 

- options available for normal 
frozen fra9ments 


Figure 1: Schematic overview of the fragment-based implementation. The 
implementation supports nonfrozen fragments, normal frozen fragments and 
frozen fragments for which the density is relaxed in freeze-and-thaw cycles. 
In addition, a number of options are available for each fragment. Reprinted 
with permission from ” A Flexible Implementation of Frozen-Density Embed¬ 
ding for Use in Multilevel Simulations”, C. R. Jacob, J. Neugebauer, and L. 
Visscher, J. Comput. Chem. 29, 1011-1018 (2008). Copyright (2015), John 
Wiley and Sons. 
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Figure 2: Schematic overview of the interlocking cells concept. The cells 
of subsystem (1) and (2) are used to evaluate the wave functions of the 
subsystems, while the supercell (E) is used to evaluate the density and the 
potential. 
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Figure 3: (a): SCF procedure used in [82] as alternative to the freeze-and- 
thaw scheme, (b): MPI architecture of a flexible and parallel Subsystem 
DFT implementation. 
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Figure 4: The 3-FDE scheme attempts to allow FDE to simulate fragments 
connected with covalent bonds, such as polipeptides. The total electron 
density is p = pi + pn — Pcap- Reprinted with permission from ’’Density- 
functional theory approach for the quantum chemical treatment of proteins”, 
C. R. Jacob and L. Visscher, J. Chem. Phys. 128, 155102 (2008). Copyright 
2015, AIP Publishing LLC. 
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(a) (b) (c) 


Figure 5: Isosurface plot of the density difference between a supermolecular 
KS calculation and a periodic subsystem DFT calculation, (a) : Methane on 
Pt(lOO), top configuration; (b): Water on Pt(lll), horizontal configuration; 
(c): Water on Pt(lll), vertical configuration. Reprinted with permission 
from ’’Periodic subsystem density-functional theory”, A. Genova, D. Ceresoli 
and M. Pavanello, J. Chem. Phys. 141, 174101 (2014). Copyright 2015, AIP 
Publishing LCC.. 
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Figure 6: OH* (H20)ii carried out with the novel FDE implementation. 
Spin-density isosurface plot (cutoff of 10“^) obtained from FDE and regular 
KS-DFT are compared. FDE produces a (correctly) localized spin-density 
on the OH* radical. 
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Figure 7: Schematic representation of the hybrid KS-Frozen environment 
method used by Hodak et al. to simulate molecular dynamics of solvated 
systems. Reprinted with permission from ’’Hybrid ab initio Kohn-Sham 
density functional theory/frozen-density orbital-free density functional the¬ 
ory simulation method suitable for biological systems”, M. Hodak, WLu, 
and J. Bernholc, J. Chem. Phys. 128, 014101 ( 2008 )]. Copyright 2015 . AIP 
Publishing LLC. 
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Figure 8: Oxygen-oxygen g{r) distribution for several combinations of 
kinetic/exchange-correlation functionals as calculated by lannuzzi et ai. 
Reprinted from ’’Density functional embedding for molecular systems”, M. 
lannuzzi, B. Kirchner, and J. Butter, Chem. Phys. Lett. 421, 16 ( 2006 ). 
Copyright ( 2015 ), with permission from Elsevier. 
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energy (eV) 

Figure 9: The electronic spectrum of the stacked benzene-fulvene dimer at 
the separation distance of 5 A, obtained using the rt-TDDFT and FDE-rt- 
TDDFT methods. 
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Figure 10: The electronic spectrum of the isolated fulvene molecule obtained 
using rt-TDDFT and the fulvene molecule in the benzene-fulvene dimer at 5 
A separation, obtained using subsystem rt-TDDFT. 
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Figure 11: The electronic spectrum of the fulvene molecule obtained using 
coupled and uncoupled subsystem rt-TDDFT methods. 
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Figure 12: Excitation energy transfer between two Na 2 molecules at the 
separation distance of 6.6 A, after the excitation of the donor molecule with 
an applied electric held with 2.18 eV frequency. 
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Approximations 

Integrand of Eq.(76) 

Acronym 

Exact 

xr'KUf 


RPA 

& Perturbative 

u.X X u.X 

Xl |r_r'|Tj 

GGP^ 

T//J 
& RPA 

& Perturbative 

Xl\r^r'\Xj 

GGP^ 

nonoverlapping 
subsystems 
& RPA 

& Perturbative 

Xl\r-r'\XJ 

GGP 


Table 3: Approximations to Eq. (76) proposed in Ref. [115] and their respec¬ 
tive acronym. “GCP” stands for Generalized Gasimir-Polder formula. 
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